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surveillance  radar  array  applications.  In  such  systems,  the  array  elements  constitute 
the  channel  outputs.  Ground  clutter.  Interference  sources,  and  noise  sources  are 
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the  technique  to  surveillance  radar  array  systems  was  established  by  identifying  several 
radar  operational  conditions  wherein  target  and/or  clutter  exhibit  non-Gaussian 
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1.0  INTRODUCTION 


Various  important  applications  involve  the  processing  of 
time-tagged  streams  of  data  collected  over  multiple  spatially- 
distributed  channels.  Such  applications  arise  in  the  military 
area  as  well  as  in  the  commercial  sector.  Military  applications 
include  radar,  sonar,  sensor  fusion,  communications,  and  seismic 
event  detection  for  monitoring  of  nuclear  test  ban  treaty 
compliance.  Non-military  applications  include  non-destructive 
evaluation,  earthquake  activity  detection,  medical  technology,  and 
flood  forecasting  and  control  in  hydrologic  systems.  The  data  can 
be  collected  either  simultaneously  over  all  channels,  or  with  a 
time  delay  between  channels.  However,  in  most  cases  the  data 
exhibits  both  temporal  and  spatial  correlations.  This 
characteristic  offers  an  opportunity  for  joint  processing  of  the 
spatial  and  temporal  domains  using  techniques  that  can  handle  such 
conditions . 

In  model-based  multichannel  detection  the  objective  is  to 
determine  the  presence  (or  absence)  of  a  desired  signal  component 
in  the  sensed  data.  This  is  carried  out  by  identifying  the 
parameters  of  an  analytic  model  fit  to  the  observed  phenomenology. 
In  this  program  the  model  class  used  to  fit  the  data  is  a  subset 
of  the  class  of  linear  systems,  and  the  techniques  used  to 
implement  the  parameter  identification  are  based  on  processing 
estimates  of  the  third-order  statistics  of  the  multichannel  data. 

The  terminology  "higher-order  statistics"  (HOS)  generally 
denotes  all  statistics  of  order  higher  than  two  (Brillinger  and 
Rosenblatt,  1967;  Mendel,  1991);  this  usage  is  adopted  in  this 
report.  Higher-order  statistics  are  of  relevance  for  data  which 
does  not  satisfy  a  Gaussian  distribution  since,  in  general, 
statistics  of  all  orders  are  required  to  describe  distributions 
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other  than  the  Gaussian  (recall  that  the  Gaussian  distribution  is 
described  completely  by  first-  and  second-order  statistics). 
Cumulants  are  a  special  type  of  statistic  with  important  practical 
implications  due  to  the  fact  that  cumulants  of  orders  third  and 
higher  are  identically  zero  for  stationary  Gaussian-dis tributed 
random  processes  (for  Gaussian-distributed  processes,  the  first- 
order  cumulant  is  the  mean  value,  and  the  second-order  cumulant  is 
the  covariance  sequence) .  Third-order  cumulants  are  of  practical 
interest  because  their  estimation  requires  significantly  less 
computations  than  the  estimation  of  cumulants  of  orders  four  and 
higher.  Besides,  in  many  cases  it  is  sufficient  to  process  only 
third-order  cumulants  in  order  to  determine  non-Gaussianity  of  the 
data,  and  to  allow  effective  estimation  of  model  parameters. 

As  indicated  above,  model-based  multichannel  detection 
methodology  and  techniques  using  higher-order  statistics  have 
applicability  in  several  areas.  However,  this  report  is  focused 
on  surveillance  radar  systems  in  general,  and  airborne 
surveillance  radar  systems  in  particular,  because  it  is  the 
application  of  most  interest  to  the  Surveillance  Technology 
Division  at  Rome  Laboratory  (RL) .  Also,  by  carrying  out  analyses 
and  simulations  in  the  context  of  a  specific  application  it  is 
possible  to  reach  a  better  understanding  of  the  utility  and 
applicability  of  the  algorithms  and  methodologies. 

The  fundamental  precondition  for  the  applicability  of  HOS  for 
model-based  radar  detection  is  the  non-Gaussian  nature  of  the 
received  radar  return,  since  HOS  vanish  in  the  case  of  Gaussian- 
distributed  processes  and  therefore  are  useless  for  model 
parameter  estimation.  An  important  outcome  of  Phase  I  is  the 
identification  of  several  surveillance  radar  operational  scenarios 
wherein  target  and/or  clutter  exhibit  non-Gaussian  statistics. 
One  such  scenario  is  the  surveillance  of  large  ocean  areas  at  low 
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grazing  angles  to  detect  long-range  cruise  missiles  flying  at  low 
altitudes  over  the  ocean.  Other  important  scenarios  are  listed  in 
Section  2.1.4. 

Utilization  of  higher-order  cumulants  in  the  context  of 
surveillance  radar  systems  is  motivated  also  by  the  fact  that 
receiver  noise  and  interference  sources  are  usually  Gaussian- 
distributed  and  statistically-independent  from  the  non-Gaussian 
clutter  and  target  return,  and  corrupt  the  radar  signal  in  an 
additive  manner.  Under  these  conditions  the  higher-order 
cumulants  of  the  signal  at  each  channel  output  are  the  cumulants 
of  the  target  and  clutter  only.  That  is,  noise  and  interference 
effects  are  eliminated  from  parameter  identification  algorithms 
based  on  higher-order  cumulants.  In  practice,  the  effects  of 
Gaussian-distributed  noise  and/or  interference  sources  are 
mitigated,  rather  than  eliminated  completely,  because  the  sequence 
used  to  estimate  the  cumulants  is  of  finite  duration. 

The  model-based  multichannel  detection  methodology  for  radar 
arrays  can  be  summarized  as  follows.  Target  detection  is 
accomplished  by  evaluating  the  output  of  two  multichannel  linear 
time- invariant  (LTI)  filters,  with  one  filter  matched  to  the 
condition  of  clutter  and  noise  only,  and  the  other  filter  matched 
to  the  condition  of  target  embedded  in  clutter  and  noise.  Each 
filter  is  a  time  series  model  fitted  to  observed  phenomenology. 
Filter  design  can  be  carried  out  either  off-line  or  on-line  in  an 
adaptive  manner.  Multichannel  filter  parameters  are  estimated 
using  algorithms  based  on  third-order  cumulants,  where  the 
cumulants  of  the  multichannel  data  are  estimated  jointly  over  the 
spatial  and  temporal  domains.  Such  algorithms  are  well  suited  for 
the  condition  of  a  desired  non-Gaussian  target  signal  embedded  in 
additive  non-Gaussian  clutter  and  additive  Gaussian  noise. 
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Other  detection  methodologies  are  suitable  for  applications 
where  non-Gaussian  signals  arise.  Of  particular  interest  is  the 
methodology  developed  by  Rangaswamy,  Weiner,  and  Michels  (1993) . 
Their  approach  addresses  innovations-based  multichannel  detection 
for  radar  surveillance  scenarios  in  which  the  spatial  behaviour  of 
the  return  signal  over  a  large  number  of  radar  resolution  (range- 
azimuth)  cells  is  described  by  a  non-Gaussian  distribution,  but 
the  temporal  behaviour  of  the  return  signal  at  each  individual 
resolution  cell  is  Gaussian.  The  spatial  return  signal  in  these 
scenarios  belongs  to  the  class  of  random  processes  known  as 
spherically  invariant  random  processes  (SIRPs).  Spatial  SIRPs 
arise  in  non-homogeneous  clutter  scenarios  where  the  temporal 
behaviour  of  the  return  signal  from  each  individual  radar 
resolution  cell  is  Gaussian-distributed,  but  the  average  return 
power  fluctuates  randomly  from  cell  to  cell.  Thus,  in  an  SIRP 
clutter  model  the  conditional  probability  density  function  (PDF) 
at  each  individual  resolution  cell  is  Gaussian. 

The  SIRP-based  detection  methodology  of  Rangaswamy,  Weiner, 
and  Michels  (1993)  is  an  optimal  formulation  for  the  case  of  a 
known,  constant  target  embedded  in  SIRP-type  clutter.  This 
methodology,  however,  cannot  be  applied  in  scenarios  where  the 
temporal  behaviour  of  the  return  signal  from  each  radar  resolution 
cell  is  described  by  a  non-Gaussian  distribution.  For  the  same 
reason,  SIRPs  cannot  be  used  to  model  the  temporal  non-Gaussian 
behaviour  of  a  target  in  a  resolution  cell . 

The  HOS-based  detection  methodology  developed  in  Phase  I  can 
be  applied  to  non-Gaussian  temporal  as  well  as  spatial  clutter 
scenarios.  Additionally,  the  HOS-based  methodology  applies  in 
cases  where  the  target  behaves  according  to  a  non-Gaussian 
distribution.  It  follows  that  the  methodology  developed  in  Phase 
I  has  broader  applicability  than  the  SIRP-based  methodology.  But 
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the  SIRP-based  methodology  is  optimal  for  the  known,  constant 
target  case.  Each  approach  has  its  merits,  and  both  deserve 
continued  consideration  and  development. 

1 . 1  Notation 


The  notation  conventions  established  herein  complement  the 
List  of  Mathematical  Symbols.  All  variables  and  functions  are  in 
Helvetica,  Chicago,  Symbol,  or  Zapf  Chancery  fonts.  Vector 
variables  are  denoted  by  underscored  lower-case  letters  (X)  , 
including  Greek  letters  (a) .  Matrices  are  denoted  by  upper-case 
letters  (A),  including  Greek  letters  (Q).  Some  scalars  are  denoted 
also  by  upper-case  letters  (for  example,  J  denotes  the  number  of 
channels  in  a  surveillance  radar  array) .  Upper-case  script 
letters  (5)  denote  sets  or  linear  spaces.  The  letter  j  is  used  to 

denote  the  imaginary  unit,  and  is  used  also  as  an  integer  index; 
the  intended  usage  is  clear  from  the  context. 

1 . 2  Report  Overview 

An  introduction  to  the  model-based  multichannel  detection 
problem  is  presented  in  Section  2.0.  This  section  also  presents  a 
discussion  of  surveillance  radar  problems  with  emphasis  on  non- 
Gaussian  operational  scenarios  (Section  2.1).  Time  series  models, 
third-order  cumulants,  and  the  associated  definitions  are 
discussed  in  Sections  2.2  and  2.3,  respectively.  The  SSC  model- 
based  multichannel  detection  methodology  is  described  in  Section 
3.0,  including  an  important  pre-processing  option  which  allows 
utilization  of  third-order  cumulants  in  cases  where  the  PDF  of  the 
radar  return  is  symmetric  or  Gaussian  (Section  3.3)  .  HOS-based 
multichannel  time  series  parameter  identification  algorithms  are 
discussed  in  Section  4.0.  Some  of  these  algorithms  have  been 
modified  and/or  extended  by  SSC  to  handle  multichannel  complex- 
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valued  sequences.  Section  5.0  discusses  results  from  simulation- 
based  analyses  carried  out  in  Phase  I  to  investigate  key  concepts 
and  issues.  This  includes  results  which  establish  the  detection 
methodology  potential  for  discrimination  between  target-present 
and  target-absent  hypotheses.  A  summary  of  the  key  technical 
issues  identified  in  Phase  I  is  presented  in  Section  6.0.  The 
main  conclusions  and  recommendations  borne  out  of  this  study  are 
discussed  also  in  Section  6.0. 


2 . 0  PROBLEM  DEFINITION 


Within  the  confines  of  airborne  surveillance  radar 
applications,  multichannel  formulations  arise  naturally  in  various 
distinct  configurations.  In  a  pre-detection  radar  fusion 
configuration  involving  two  or  more  radar  systems  irradiating  the 
same  physical  area,  each  channel  is  the  output  of  an  individual 
radar  with  a  particular  set  of  operational  parameters.  The 
radiated  frequency  and/or  any  other  operational  parameters  (such 
as  polarization)  may  be  different  from  radar  to  radar  (channel  to 
channel) .  In  a  radar  array  configuration  each  channel  is  the 
output  of  an  individual  antenna  element,  or  the  output  of  an 
individual  subarray  (a  subarray  is  a  collection  of  antenna 
elements) .  The  methodologies  presented  herein  can  be  applied  to 
these  two  configurations  as  well  as  to  others.  Notwithstanding, 
the  radar  array  configuration  is  selected  as  the  baseline  for 
discussions . 

For  the  radar  applications  of  interest  herein,  radar  target 
detection  can  be  posed  as  a  hypothesis  testing  problem  involving 
the  following  two  hypotheses: 

Hq :  Desired  target  signal  is  absent 

H1 :  Desired  target  signal  is  present 

where  H0  is  the  null  hypothesis,  and  Hi  is  the  al  t  ernat  ive 
hypothesis .  The  radar  detection  processor  must  select  between 
these  two  hypotheses.  Selection  of  Hi  when  Hi  is  true  is  a  correct 
detection,  and  a  common  measure  of  radar  detection  performance  is 
the  probability  of  detection.  Two  types  of  detection  errors  can 
occur:  (1)  selection  of  Hi  when  H0  is  true,  which  results  in  a 
false  alarm,  and  (2)  selection  of  H0  when  Hi  is  true,  which  results 
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in  a  missed  detection.  Each  of  these  two  errors  has  an  associated 
probability  measure:  probability  of  false  alarm  for  a  type  1 
error,  and  probability  of  miss  for  a  type  2  error. 

Model-based  multichannel  detection  techniques  of  the  type 
pursued  herein  are  based  on  the  premise  that  the  received  radar 
time  sequence  (after  demodulation  to  baseband  and  temporal 
sampling)  can  be  represented  accurately  for  each  hypothesis  as  the 
output  of  an  analytic  model  defined  by  a  unique  set  of  model 
parameters.  It  is  assumed  also  that  the  two  parameteric 

representations  (one  corresponding  to  each  hypothesis)  are 
sufficiently  distinct  to  allow  selection  of  the  correct  hypothesis 
by  evaluation  of  measures  that  are  sensitive  to  those  differences. 
Linear  time-invariant  (LTI)  models  are  adopted  in  this  program 
because  such  models  allow  significant  modeling  freedom  while 
presenting  manageable  analytical  and  computational  burdens.  An 
important  class  of  LTI  models  is  the  so-called  time  series  model 
class,  which  includes  moving-average  (MA)  ,  auto-regressive  (AR) , 
and  auto-regressive  moving-average  (ARMA)  models;  this  is  the 
model  class  adopted  in  this  program.  The  model  output  is 

corrupted  by  additive  Gaussian-distributed  noise,  which  is  assumed 
to  be  uncorrelated  in  time  (white) . 

In  the  philosophy  summarized  above,  a  time  series  model  is  a 
representation  model  for  the  spatial  or  temporal  statistics  of  the 
received  radar  signal;  that  is,  a  model  which  fits  experimental 
and/or  simulated  statistical  data.  In  contrast,  a  physical  model 
is  the  result  of  a  theoretical  formulation  of  the  physical 
phenomenology.  Of  course,  physical  models  also  fit  real  data. 
Haykin,  Currie,  and  Kesler  (1982),  and  Metford  and  Haykin  (1985) 
have  demonstrated  that  time  series  models  provide  an  excellent  fit 
to  surveillance  radar  data.  Their  results  motivate,  in  part,  the 
research  in  this  program. 


A  particular  measure  that  has  produced  robust  analytical  and 
experimental  results  in  the  model-based  detection  context  is  the 
log-likelihood  ratio  (LLR)  test.  This  test  is  the  outcome  of 
solving  the  hypothesis  testing  problem  using  the  Neyman-Pearson 
criterion.  The  Neyman-Pearson  criterion  maximizes  the  probability 
of  detection  for  a  given  level  of  probability  of  false  alarm.  The 
LLR  test  in  the  context  of  model-based  detection  is  calculated 
using  the  innovations  sequence  at  the  output  of  each  of  two  LTI 
filters,  one  for  each  hypothesis  (of  course,  only  the  filter 
output  corresponding  to  the  true  hypothesis  is  a  true  innovations 
sequence).  Metford  and  others  (Metford  et  al .  ,  1982;  Metford 
(1984);  Metford  and  Haykin,  1985)  have  formulated  the  discrete¬ 
time  version  of  the  LLR  detection  test  using  model-based 
innovations,  and  labeled  the  methodology  as  the  innovations  based 
detection  algorithm  (IBDA) . 

In  many  surveillance  scenarios  the  received  radar  time  series 
includes  at  least  one  component  (either  target  or  clutter) 
described  by  a  non-Gaussian  distribution,  and  corrupted  with 
additive  Gaussian-dis tributed  receiver  and/or  interference  noise. 
An  important  feature  which  distinguishes  the  Gaussian  distribution 
from  others  is  the  fact  that  the  Gaussian  distribution  is  defined 
completely  in  terms  of  its  first-  and  second-order  statistics. 
Thus,  statistics  of  order  third  and  higher  assume  an  important 
role  when  non-Gaussian  distributions  are  involved.  The  cumulants 
of  a  distribution  is  a  set  of  statistics  of  particular  interest 
because  higher-order  cumulants  are  zero  for  the  Gaussian 
distribution,  but  non-zero  for  non-Gaussian  distributions  except 
that  odd-ordered  cumulants  are  zero  for  distributions  with 
symmetry  about  the  mean  value  (see,  for  example,  Mendel  [1991]) . 
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The  above  discussion  suggests  that  parameter  identification 
techniques  based  on  cumulants  can  be  used  for  model-based 
multichannel  detection  in  non-Gaussian  scenarios.  In  the  approach 
pursued  herein  the  parameters  of  a  LTI  time  series  model  (MA,  AR, 
or  ARMA)  are  identified  using  techniques  based  on  processing 
estimates  of  third-order  cumulants.  The  effects  of  additive 
Gaussian  noise  are  mitigated  since  the  third-order  cumulants  of  a 
Gaussian-distributed  process  are  zero.  Third-order  cumulants  are 
preferred  because  their  estimation  from  data  requires  fewer 
computations  than  the  estimation  of  cumulants  of  order  four  and 
higher.  Parameter  identification  algorithms  are  developed  based 
on  the  work  of  Giannakis,  Inouye,  and  Mendel  (1989)  and  Raghuveer 
(1986) ,  modified  and  extended  appropriately  to  handle  the  case  of 
complex-valued  data.  An  optimal  LLR  detection  rule  for  the 
outputs  of  the  LTI  filters  associated  with  each  hypothesis 
requires  knowledge  of  the  multivariate  PDF  under  each  of  the 
hypotheses.  In  most  non-Gaussian  cases  this  information  is 
unavailable.  However,  Metford  (1984)  has  demonstrated  that  the 
LLR  detection  rule  for  the  Gaussian  case  is  a  suboptimal, 
approximate  detection  rule  for  the  non-Gaussian  case.  This 
approximation  was  adopted  for  the  work  reported  herein. 

Strictly  speaking,  the  utility  of  third-order  cumulants  is 
restricted  to  cases  where  the  PDF  of  the  radar  return  (quadrature 
components)  is  asymmetric  about  the  mean  value.  However,  it 
appears  that  the  cases  where  the  PDF  is  symmetric  can  be  addressed 
by  first  applying  a  zero-memory  nonlinear  transformation  to  the 
data.  Such  an  approach  has  been  proposed  recently  by  Zheng, 
McLaughlin,  and  Mulgrew  (1993),  with  the  logarithm  function  as  the 
zero-memory  nonlinear  transformation.  Their  approach  is  discussed 
in  Section  3.3.  Results  from  simulation-based  analyses  carried 
out  in  Phase  I  indicate  the  logarithm  is  indeed  useful  for 
introducing  asymmetry  in  the  data  PDF  (see  Section  5.1)  .  This 
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issue  has  practical  (as  well  as  theoretical)  significance  because 
in  many  radar  systems  the  logarithm  function  is  used  to  reduce,  the 
dynamic  range  of  the  radar  return. 

2 . 1  Airborne  Surveillance  Radar 

Consider  a  coherent  radar  array  with  J  spatial  channels  (each 
channel  is  the  output  of  either  an  individual  array  element  or  a 
sub-array  composed  of  multiple  array  elements),  as  indicated  in 
Figure  2-1.  The  radar  is  coherent  because  both  in-phase  and  out- 
of-phase  quadrature  components  are  generated  in  each  channel.  In 
a  surveillance  scenario  involving  a  down-looking  radar  onboard  an 
aircraft  (see,  for  example,  Jaffer  et  al . ,  [1991],  or  Rangaswamy 
et  al .  ,  [1993]),  the  J-element,  discrete- time ,  baseband, 
stationary,  complex-valued,  finite-duration,  vector  sequence  {xk(n)  I 
n  =  0,  1,  .  .  .  ,  N-1}  in  Figure  2-1  is  the  return  from  the  kth  radar 
resolution  (range- azimuth)  cell  received  at  the  J  channels  for  the 
duration  of  the  coherent  processing  interval  (CPI) . 


{xk(n)} 


xk(n)  :  complex-valued,  J-dimensional  vector  of  the  return  from 
the  kth  range  gate  corresponding  to  the  nth  pulse 


Figure  2-1.  Multichannel  signal  in  a  coherent  surveillance  radar 

array . 


In  the  context  of  a  hypothesis  testing  formulation  applied  to 
the  kth  resolution  cell,  the  null  hypothesis  corresponds  to  the 
case  of  target  absent;  the  alternative  hypothesis  corresponds  to 
the  case  of  target  present.  Under  the  null  hypothesis,  the 
baseband  vector  process  tek(n)}  contains  clutter,  interference,  and 
noise.  Under  the  alternative  hypothesis,  {xk(n)}  also  contains 
target  information.  This  can  be  summarized  as 


(2-la) 

H0: 

><k(n) 

=  ck(n)  +  wk(n) 

0  <n  <  N-1 

(2-lb) 

Hi. 

Xk(n) 

=  sk(n)  +  ck(n)  +  wk(n) 

0  <  n  <  N-1 

where  {ck(n)}  denotes  the  clutter  process,  {wk(n)}  denotes  the 

combination  of  all  noise  (receiver  noise,  etc.)  and  interference 
processes,  and  {sk(n)}  denotes  the  desired  target  signal  process. 
Table  2-1  lists  the  conditions  assumed  on  the  vector  process  {xk(n)}. 
Notice  that  selection  of  the  initial  time  as  n  =  0  in  Equations  (2- 
1)  implies  no  loss  of  generality  since  the  process  is  assumed  to 
be  stationary  (in  practice  it  suffices  that  the  process  be 
stationary  at  least  for  the  duration  of  the  CPI) . 

In  a  typical  surveillance  application  over  a  large  number  of 
resolution  cells,  the  data  from  each  radar  resolution  cell  is 
processed  and  a  decision  is  made  as  to  the  presence  or  absence  of 
a  target.  Various  resolution  cells  can  be  processed  in  parallel 
or  sequentially,  depending  on  the  design  parameters  of  the  radar 
system.  Consider  the  data  {xk(n)}  from  one  CPI  for  a  specific 

resolution  cell  to  which  the  detection  criterion  is  being  applied. 
This  data  set  is  referred  to  as  the  primary  data .  Returns  from  G 
cells  in  the  neighborhood  of  the  kth  cell  are  described  by  the  same 
probability  distribution  as  the  kth  cell,  but  it  is  assumed  that  a 
target  is  not  present  in  those  cells.  It  is  further  assumed  that 
the  data  in  those  cells  is  statistically  independent  of  the  data 
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in  the  kth  cell.  The  data  from  these  G  range  cells  is  referred  to 
as  the  secondary  data,  and  is  used  to  set  the  detection  threshold 
for  the  kth  resolution  cell.  As  the  detection  process  continues, 
the  current  primary  data  becomes  part  of  the  secondary  data 
corresponding  to  a  different  resolution  cell. 


BRSEBRND  UECTOR  SEQUENCE  ASSUMPTIONS  UNDER  BOTH  HYPOTHESES 

•  Zero  Mean 

•  Stationary 

•  Ergodic 

•  Nyquist  Sampling  Criterion  Is  Satisfied 

•  Noise  Process  Is  Gaussian-Distributed  and  White 

•  Target  Process  Is  Non-Gaussian 

•  Clutter  Process  Is  Non-Gaussian 

•  Target,  Noise,  and  Clutter  Processes  fire  Independent 

•  Can  Be  Modeled  as  Output  of  a  LTI  System 


Table  2-1.  Baseband  vector  sequence  assumptions  under  both 

hypotheses . 

From  this  point  on,  the  subscript  k  on  the  baseband  vector 
sequence  (this  subscript  denotes  the  resolution  cell)  will  be 
deleted  to  simplify  notation.  This  is  appropriate  also  because, 
as  discussed  above,  the  data  from  each  resolution  cell  assumes  the 
status  of  primary  data  sometime  during  the  process . 

In  most  of  the  radar  technical  literature  the  terminology 
"Gaussian-distributed  signal"  denotes  a  coherent  radar  return  with 
Gaussian-distributed  quadrature  components,  Rayleigh-distributed 
envelope,  and  exponentially-distributed  power  (or  radar  cross- 
section,  RCS) .  Correspondingly,  "non-Gaussian  signal"  refers  to  a 
radar  return  with  either  non-Rayleigh  envelope  or  non-exponential 
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power.  This  terminology,  although  common,  can  be  misleading  in 
the  context  of  algorithms  and  methodologies  based  on  third-order 
cumulants  because  the  shape  of  the  distribution  of  the  quadrature 
components  is  a  key  issue  when  third-order  cumulants  are  involved 
(as  well  as  all  odd-ordered  cumulants  of  higher  orders)  .  To 
illustrate  this  point  consider  the  case  of  a  signal  consisting  of 
a  constant-frequency  sinusoid  (representing  a  constant-velocity 
moving  target)  in  additive  Gaussian  white  noise.  In  this  case  the 
combined  signal  (target  plus  noise)  envelope  is  described  by  the 
Rice  distribution,  the  phase  is  not  uniformly-distributed,  and  the 
envelope  and  phase  are  statis tically-dependent  (Thomas,  1969). 
This  case  is  referred  to  in  the  literature  as  non-Gaussian,  but 
such  a  statement  is  misleading  because  the  quadrature  components 
are  jointly  Gaussian-distributed  and  independent  of  each  other, 
and  the  marginal  distribution  of  each  is  Gaussian  with  non-zero 
mean.  Estimation  of  the  higher-order  cumulants  of  the  quadrature 
components  is  carried  out  after  subtracting  the  mean,  which 
results  in  zero-mean  Gaussian-distributed  quadrature  components. 
Thus,  for  a  Rice-distributed  envelope  all  higher-order  cumulants 
of  the  quadrature  components  are  zero  (after  mean  removal) . 

2.1.1  Alternative  Probabilistic  Descriptions 

A  complete  probabilistic  model  for  radar  data  consists  of  the 
multivariate  PDF  of  the  received  radar  time  sequence  for  the  full 

CPI,  {x(n)  I  n  =  0,  1 . N} .  The  multivariate  PDF  model  inherently 

includes  the  temporal  and  spatial  statistical  information 
(moments,  etc.)  required  to  estimate  the  parameters  of  a  time 
series  model  to  represent  the  data.  A  multivariate  PDF  model  is 
required  also  under  each  of  the  two  hypotheses  in  order  to 
determine  the  optimal  LLR  detection  rule.  Unfortunately,  a 
temporal  multivariate  PDF  model  is  available  only  for  a  limited 
number  of  scenarios  and  conditions  that  result  in  Gaussian- 
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distributed  return  signals .  In  most  cases  involving  non-Gaussian- 
distributed  return  signals  the  available  information  is  a 
univariate  PDF  model  for  either  the  envelope  or  the  power  of  each 
scalar  element  of  {x(n)}  at  each  discrete  instant  of  time.  Some  of 
the  available  univariate  PDF  models  are  physical  models,  whereas 
others  are  representation  models . 

Multivariate  as  well  as  univariate  PDF  descriptions  can  be 
represented  in  either  one  of  two  ways:  as  a  function  of  complex¬ 
valued  quantities,  or  as  a  function  of  real-valued  quantities. 
Issues  associated  with  the  two  representations  are  considered  in 
this  section,  with  emphasis  on  those  issues  that  have  an  impact  on 
model-based  multichannel  detection  for  non-Gaussian  cases. 

Let  the  complex-valued  baseband  process  {x(n)}  be  defined  in 
terms  of  its  real  and  imaginary  components  as 

(2-2)  x(n)  =  x,(n)  +  j  Xq(ii) 

where  the  subscripts  I  and  Q  denote  in-phase  and  out-of -phase , 
respectively.  The  real-valued  vector  sequences  {Xj(n)}  and  {xQ(n)}  are 

the  quadrature  components  of  {x(n)}  •  Since  the  process  {x(n)}  is 
zero-mean,  both  of  these  vectors  are  zero-mean  also.  The  auto¬ 
covariance  and  cross-covariance  sequences  of  these  two  components 
are  defined  as 

(2-3)  RM(m)  =  E[x,(n)x[(n-m)] 

(2-4)  RQQ(m)  =  E[xQ(n)x^(n-m)] 

(2-5)  RIQ(m)  =  E[  x,(n)  x^(n-m)] 

(2-6)  RQI(m)  =  E[xQ(n)x[(n-m)] 
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Auto-covariances  R|i(m)  and  RQo(m)  are  both  symmetric  and  (at  least) 
positive  semi-definite.  Cross-covariances  Rio(m)  and  RQ|(m)  do  not 
exhibit  any  particular  structure,  but  it  is  obvious  from  Equations 
(2-5)  and  (2-6)  that 

(2-7)  R|Q(m)  =  RQ,(m) 


The  covariance  sequence  of  the  process  {x(n)}  is  defined  as 
(2-8a)  RXx(m)  =  E  [x(n)  xH(n-m)  ] 

(2-8b)  Rxx(m)  =  R||(m)  +  RQQ(in)+jRQ|(m)-jR|Q(m) 

Matrix  Rxx(m)  is  Hermitian  and  (at  least)  positive  semi-definite. 
Also,  the  following  relationship  is  true 

(2-9)  Rxx(-m)  =  RxX(m)  Vm 


This  covariance  sequence  constitutes  the  second-order  statistics 
for  the  process  {x(n)} .  Since  the  process  has  mean  zero,  the  matrix 
elements  of  the  covariance  sequence  are  the  second-order  cumulants 
of  the  process . 


Now  define  three  JN  -dimensional  vectors  by  concatenating  all 
the  J-dimensional  vectors  in  the  complete  CPI  as  follows: 


(2-10) 


x,(0;N-1) 


X|(0) 

x,(1) 


Lx,(N-1)J 
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"  xQ(0) 

XqO) 

(2-11)  xq(0;N-1)  =  . 

_xq(N-1) 

x(O) 
x(1) 

(2-12)  x(0;N-1)  =  x,(0;N-1)  +  jxQ(0;N-1)  = 

.  x(N-1) J 

The  vectors  in  Equations  (2-10 )-( 2-12 )  contain  all  the  data  for  a 
specific  resolution  cell  for  one  CPI.  Define  also  a  2JN- 
dimensional,  real-valued  vector  as 

x,(0;N-l) 

(2-13)  r(0;N-1)  = 

xq(0;N-1) 

Since  the  process  {x(n)}  is  zero-mean,  all  of  the  vectors  in 
Equations  (2-10)  -  (2-13 )  are  zero-mean.  Vector  x(0;N-1)  contains  the 
same  information  as  vector  r(0;N-1);  however,  it  turns  out  that  the 
probabilistic  description  of  the  baseband  process  based  on  the 
complex- valued  representation  of  the  information  is  more 
restrictive . 

First  it  is  convenient  to  introduce  additional  definitions. 
The  covariance  matrix  of  the  JN-dimensional  vector  x(0;N-1),  is 
defined  as 

% 

'  Rxx(°)  RxxH)  RxxO-N) 

(2-14)  ^NN  =  E[x(0;N-1)xH(0;N-1)]  =  Rx*(1)  Rxx(0)  Rxx(2'N) 

_Rxx(N-D  Rxx(N-2)  •••  Rxx(0) 
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Matrix  N  is  Hermitian  and  block  Hermitian,  and  it  exhibits  block 

Toeplitz  structure  (the  ( i , j ) th  block  element  of  a  block  Toeplitz 
matrix  is  a  function  of  i-j)  .  The  block  elements  of  N  are 

matrix  elements  of  the  covariance  sequence  {Rxx(m)} .  Also  relevant 
is  the  covariance  matrix  of  r(0;N-1),  which  is  defined  as 


(2-15a)  Rnn  =  E[  r(0;N-1)  rT(0;N-1)  ] 


R„(0)  ■ 

■■  R„(1-N) 

R|q(0) 

■■  R,q(1-N) 

R„(N-1)  • 

R„(0) 

R,q(N-1)  • 

R|q(0) 

RQ|(0) 

•  RQ|(1-N) 

a  ■  ■ 
o 

DC 

■  Rqq(i-n) 

Rqi(N-I)  • 

■  •  RQ|(0) 

Rqq(N-I)  .. 

-  rqq(0) 

Matrix  RN N  is  symmetric  and  block  symmetric;  however,  in  general 

it  does  not  exhibit  block  Toeplitz  structure.  The  block  elements 
of  Rn  n  are  the  real  and  imaginary  matrix  components  of  the 

covariance  sequence,  Equation  (2-8b) . 

Having  established  the  relevant  notation  and  correlation 
matrix  definitions,  the  concept  of  circular  symmetry  can  be 
defined.  A  complex-valued  vector  process  {x(n)}  with  covariance 
sequence  defined  as  in  Equation  (2-8)  is  circularly  symmetric  if 
the  following  conditions  are  satisfied  (Goodman,  1984;  Papoulis, 
1984)  : 


(2-16) 

E[x,(n)]  =  E[XQ(n)]  =  0 

V  n 

(2-17a) 

Rn(m)  =  RQa(m) 

V  m 
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(2-i7b)  R|Q(m)  =  -  RQ|(m) 


V  m 

Condition  (2-16)  requires  the  radar  return  to  have  mean  zero, 
which  is  satisfied  by  assumption  (see  Table  2-1) ;  Conditions  (2- 
17)  place  structural  constraints  on  the  covariance  sequence  of 
{x(n)} .  Notice  that  circular  symmetry  involves  only  first-  and 
second-order  statistics.  Thus,  most  results  available  in  the 
literature  on  circular  symmetry  involve  the  Gaussian  distribution 
(recall  that  the  Gaussian  distribution  is  defined  completely  by 
first-  and  second-order  statistics) .  One  example  of  such  results 
is  the  following:  if  a  Gaussian-distributed  baseband  sequence  {x(n)} 
is  circularly  symmetric,  then  its  associated  bandpass  process  {x(t)} 
(see  Figure  2-1)  is  a  stationary  process  (Papoulis,  1984).  ;  The 
converse  statement  is  also  true.  An  equivalent  result  for  non- 
Gaussian  processes  is  unavailable;  it  seems  that  such  a  result  may 
not  be  possible  for  non-Gaussian  processes  because  of  the  role 
played  by  higher-order  statistics  in  such  processes. 

The  probabilistic  description  of  the  complex-valued  vector 
process  {x(n)}  can  be  represented  using  either  one  of  two  analytic 
forms.  In  the  first  form  the  multivariate  PDF  is  expressed  in 
terms  of  the  JN-dimensional  complex-valued  vector  x(0;N-1),  with 
corresponding  complex-valued  moments.  Alternatively,  in  the 
second  form  the  multivariate  PDF  is  expressed  in  terms  of  the  2JN- 
dimensional  real-valued  vector  l(0;N-1),  with  corresponding  real¬ 
valued  moments.  It  is  well  known  that  when  {x(n)}  is  Gaussian- 
distributed,  the  two  PDF  representations  are  equivalent  if  the 
process  {x(n)}  is  circularly  symmetric  (see,  for  example,  Goodman 
[1963]) .  The  above  comments  and  conclusions  are  valid  also  for 
univariate  PDF  models  for  each  element  of  the  random  vector  x(n)  at 
each  time  instant  n. 
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Circular  symmetry  is  not  satisfied  in  many  cases,  even  when 
the  quadrature  components  are  Gaussian-distributed .  As  an  example 
consider  a  process  represented  by  the  Rice  envelope  PDF  model. 
The  quadrature  components  corresponding  to  a  Rice  envelope  process 
are  Gaussian-distributed,  but  the  mean  is  non-zero.  Since  the 
process  fails  Condition  (2-16),  it  is  not  circularly  symmetric. 
Consequently,  this  process  cannot  be  represented  with  a  PDF  model 
based  on  the  complex-valued  formulation.  Similar  behaviour  can 
occur  for  non-Gaussian  processes.  Furthermore,  even  if  a  non- 
Gaussian  process  is  circularly  symmetric,  there  may  not  be  a  valid 
complex- valued  PDF  model.  In  summary,  the  multivariate  PDF 
representation  based  on  the  real-valued  vector  r(0;N-1)  and  its 
moments  is  more  general  than  the  representation  based  on  the 
complex-valued  vector  x(0;N-1)  and  its  moments  (Williams  and 
Johnson,  1993;  Michels,  1991).  This  fact  is  significant  because 
optimal  LLR  detection  algorithms  are  developed  based  on  the 
multivariate  PDF  model  for  the  sequence. 

Based  on  the  discussion  presented  in  this  section,  it  is 
reasonable  to  adopt  the  real-valued  probabilistic  representation 
in  the  formulation  of  model  parameter  identification  algorithms 
and  detection  methodologies  for  scenarios  involving  non-Gaussian 
radar  returns.  Nevertheless,  the  methodology  and  equations  were 
developed  in  Phase  I  using  the  complex-valued  probabilistic 
representation  because  the  form  of  the  resulting  equations  and 
diagrams  are  applicable  to  both  cases  (in  general,  a  complex- 
formulation  equation  reduces  to  its  real-valued  counterpart  by 
setting  the  imaginary  components  to  zero  and  replacing  the 
Hermitian  operator  with  the  transpose  operator) .  Additionally, 
the  notation  for  the  complex-valued  case  is  more  compact  than  for 
the  real-valued  case  (compare,  for  example,  Equation  (2-14)  with 
Equation  (2-15)). 


20 


2.1.2 


Non-Gaussian  Clutter  Scenarios  and  Probabilistic  Models 


Airborne  surveillance  radars  generally  search  for  moving 
targets  over  large  areas  of  ground  surface  consisting  of  land,  or 
water  (sea,  lake,  etc.),  or  both  (coastal  regions,  etc.).  Thus, 
the  received  radar  signal  contains  land  and/or  water  clutter,  and 
the  probabilistic  description  of  these  two  clutter  types  is  of 
relevance.  The  methodology  developed  in  this  program  is  best 
suited  for  cases  where  the  return  signal  is  non-Gaussian  (unless  a 
transformation  is  applied  to  the  data  first);  thus,  identification 
of  radar  scenarios  and/or  conditions  which  result  in  non-Gaussian 
clutter  is  a  fundamental  issue. 

As  mentioned  in  Section  2.1.1,  the  multivariate  PDF  of  the 
received  radar  time  sequence  for  one  full  CPI  is  a  complete 
probabilistic  model  for  radar  data.  However,  in  most  cases  of 
non-Gaussian  clutter  the  available  PDF  model  is  a  univariate  PDF 
model  for  either  the  envelope  or  the  power  (equivalently,  radar 
cross-section)  at  the  output  of  each  individual  channel  and  at 
each  discrete  instant  of  time.  The  univariate  PDF  is  complemented 
by  temporal  correlation  information.  During  Phase  I  of  this 
program  a  literature  review  was  carried  out  to  identify  the 
various  non-Gaussian  clutter  univariate  PDF  models  and  the 
conditions  and/or  scenarios  where  such  non-Gaussian  clutter 
arises .  A  summary  of  the  most  common  univariate  PDF  models  is 
presented  in  Table  2-2.  The  last  column  in  Table  2-2  describes 
the  clutter  type  and  the  conditions  under  which  the  PDF  model  is 
valid.  Several  relevant  references  are  provided  also  in  the  last 
column.  Most  of  the  references  given  in  this  column  are  texts 
which  discuss  clutter  issues  in  detail  and  provide  extensive  lists 
of  references  for  the  original  experimental  and/or  analytical  work 
on  which  the  results  summarized  herein  are  based. 
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UNIVARIATE 
PDF  MODEL 

MODEL 

PARAMETER 

VARIATION 

TYPE 

CLUTTER  TYPE 
DESCRIPTION 

Rice 

(Non-central 

gamma) 

Envelope 

(RCS) 

Temporal 

Terrain  with  large,  fixed  scatterers  (such  as 
man-made  structures)  and  a  large  number 
of  small,  moving  scatterers  (Nathanson, 
1991;  Simkins,  1981). 

K-distribution 

Envelope 

Spatial 

Sea  clutter  at  microwave  frequencies 
(Jakeman  and  Pusey,  1976;  Nathanson, 
1991). 

K-distribution 

Envelope 

Spatial 

Polarimetric  scattering  from  certain  types  of 
terrain  (Raghavan,  1991). 

Gamma 

RCS 

Spatial 

Sea  clutter  at  S-band  and  C-band  radar 
frequencies;  used  in  RADC  model 
(Simkins,  1991). 

Log-normal 

RCS 

Temporal;  Spatial 

Sea  clutter  for  high-resolution  radar  (short- 
duration  pulses  and  narrow  beamwidth), 
low  grazing  angles  ( <  5  deg),  and  high  sea 
states  (Nathanson,  1991;  Skolnik,  1980). 

Log-normal 

RCS 

Spatial 

Various  types  of  terrain  clutter  at  low 
grazing  angles  (<  3  deg);  also  non- 
homogeneous  (composite)  terrain  (Currie, 
1989;  Nathanson,  1991;  Skolnik,  1980). 

Log-normal 

Envelope 

Temporal 

Rain  clutter  with  millimeter  wave  radar  at 
rain  rates  in  the  1  mm/hr  to  60  mm/hr  range 
(Skolnik,  1980). 

Weibull 

Envelope 

Temporal;  Spatial 

Sea  clutter  for  high-resolution  radar  (short- 
duration  pulses  and  narrow  beamwidth), 
low  grazing  angles  (  <  5  deg),  and  high  sea 
states  (Nathanson,  1991;  Skolnik,  1980). 

Weibull 

Envelope 

Spatial 

Various  types  of  terrain  clutter  at  low 
grazing  angles  (<  3  deg),  also  non- 
homogeneous  (composite)  terrain  (Currie, 
1989;  Nathanson,  1991;  Skolnik,  1980;  ). 

Weibull 

Envelope 

Spatial 

Weather  clutter  under  stormy  conditions  at 
L-band  (1.3  GHz)  over  Tokyo,  Japan 
(Sekine  et  al.,  1979). 

Table  2-2.  Non-Gaussian  clutter  fluctuation  PDF  models. 

The  univariate  PDF  models  in  Table  2-2  provide  only  part  of 
the  information  required  for  the  development  of  optimal  detection 
methodologies  and  analyses.  Minimum  additional  model  information 
required  for  suboptimal  techniques  is  a  description  of  the 
temporal  correlation  structure.  The  most  commonly  adopted 
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temporal  correlation  models  are  those  which  bound  temporal 
correlation  behaviour  as  defined  by  Marcum  and  Swerling  (1960) : 
(a)  pulse- to-pulse  fluctuation,  and  (b)  scan-to-scan  (or  pulse 
train- to-pulse  train)  fluctuation.  In  multichannel  systems 
spatial  correlation  models  are  required  also.  Analogous  to  the 
temporal  case,  the  two  models  that  bound  spatial  fluctuation 
behaviour  are:  (a)  channel- to-channel  independence,  and  (b) 
channel- to-channel  full  correlation.  Multichannel  systems  offer 
more  degrees  of  freedom,  with  the  associated  performance 
improvement  potential  and  increased  analytical  complexity. 

A  clutter  PDF  model  represents  either  a  temporal  or  a  spatial 
variation  for  the  model  parameter  (envelope  or  RCS) ,  as  indicated 
in  Table  2-2.  A  temporal  PDF  model  describes  the  instantaneous 
fluctuations  of  the  envelope  (or  power)  of  the  baseband  radar 
signal  {x(n)}  at  each  individual  resolution  cell.  A  spatial  PDF 
model  describes  cell-to-cell  fluctuations  of  the  mean  envelope  or 
mean  power  (equivalently,  RCS)  over  a  large  area  covering  many 
resolution  cells,  and  these  fluctuations  behave  according  to  the 
listed  distribution,  except  in  the  case  of  the  K-dis tribution . 
The  cell-to-cell  fluctuations  in  the  K-dis tribution  model  behave 
according  to  the  (un-normalized)  chi  distribution.  In  the  model 
referred  to  in  the  table,  the  K-dis tribution  describes  the  joint 
fluctuation  behaviour  of  the  radar  return  from  all  the  resolution 
cells  over  a  CPI  (this  is  discussed  in  more  detail  below) .  Thus, 
if  the  same  naming  convention  were  to  be  used  in  all  cases,  the 
correct  name  for  this  model  would  be  "chi".  However,  the  K- 
distribution  name  for  the  model  is  entrenched. 

In  most  spatial  model  cases  the  instantaneous  envelope 
fluctuations  in  each  individual  radar  resolution  cell  are 
Rayleigh-distributed,  and  the  instantaneous  power  fluctuations  are 
exponentially-distributed.  This  is  the  Gaussian  case.  Thus,  most 
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spatial  PDF  models  represent  non-homogeneous  clutter  with 
Gaussian-distributed  behaviour  at  each  individual  resolution  cell. 
The  spatial  vs.  temporal  distinction  is  important  in  the  context 
of  higher-order  cumulants  since  all  higher-order  cumulants  are 
zero  for  Gaussian-distributed  clutter  returns.  An  HOS-based 
detection  formulation  can  be  defined  for  each  case  (spatial; 
temporal),  as  discussed  in  Section  3.2. 

With  exception  of  the  Rice  envelope  model,  each  of  the  PDF 
models  listed  in  Table  2-2  is  incomplete  in  the  sense  that 
information  is  unavailable  in  the  literature  which  is  essential  to 
determine  the  unqualified  feasibility  of  time  series  model 
parameter  estimation  using  third-order  cumulants  of  the  quadrature 
components  ("unqualified"  implies  without  the  pre-processing  step 
introduced  in  Section  3.3) .  Specifically,  the  following 
information  is  unavailable  for  each  of  the  listed  PDF  models 
(except  Rice  envelope) : 

(a)  joint  distribution  of  the  envelope  and  phase, 

(b)  marginal  distribution  of  the  phase, 

(c)  joint  distribution  of  the  quadrature  components,  and 

(d)  marginal  distribution  of  the  quadrature  components . 

Of  these  four  items,  the  key  model  information  (in  the  context  of 
third-order  cumulants  of  the  quadrature  components)  is  the 
marginal  distribution  of  the  quadrature  components.  Given  the 
marginal  PDF,  the  asymmetry  of  the  distribution  about  the  mean 
value  can  be  determined.  An  effective  measure  of  univariate  PDF 
asymmetry  is  the  coefficient  of  skewness.  r|3,  which  is  defined  as 

(2-18)  n3  = 

CT3 

where  JI3  is  the  third  moment  about  the  mean  and  a  is  the  standard 

deviation  of  the  PDF  model  (Hastings  and  Peacock,  1974)  .  For  a 
symmetric  distribution,  pi 3  =  0 ;  a  highly  asymmetric  PDF  has  a  large 
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coefficient  of  skewness.  Variants  of  the  definition  in  Equation 
(2-18)  have  been  introduced  in  the  literature  for  specific 
objectives.  Maffett  (1989)  uses  a  modified  coefficient  of 
skewness,  where  the  modification  consists  of  replacing  a3  in  the 
denominator  by  the  product  |1G2,  with  |i  the  mean  of  the  process. 

Often  it  is  assumed  that  the  clutter  return  phase 
corresponding  to  the  temporal  distributions  in  Table  2-2  is 
uniformly-distributed  and  independent  of  the  envelope  (or  RCS)  . 
Sometimes  these  two  assumptions  are  made  to  simplify  the 
generation  of  simulated  data.  Other  times  they  are  introduced  in 
order  to  simplify  analyses.  However,  the  validity  of  these  two 
assumptions  has  not  been  justified  via  experimental  and/or 
analytical  means.  In  fact,  the  Rice  envelope  is  an  example  of  the 
opposite:  in  the  Rician  envelope  model  the  phase  and  the  envelope 
are  dependent,  and  the  phase  marginal  PDF  is  not  uniform  (Thomas, 
1969) .  The  quadrature  components  of  the  Rice  envelope  model  are 
independent  and  Gaussian-distributed,  though.  The  validity  (or 
lack  thereof)  of  the  two  assumptions  on  phase  (uniformly- 
distributed  and  independence  of  phase  and  envelope)  for  any  of  the 
non-Gaussian  cases  in  Table  2-2  (except  Rice  envelope)  has 
significant  implications  in  the  context  of  third-order  cumulants . 
This  is  so  because  it  can  be  shown  that  if  the  phase  is  uniformly- 
distributed  and  independent  of  the  envelope  (or  RCS)  ,  then  the 
marginal  PDF  of  each  of  the  two  quadrature  components  is  symmetric 
about  the  mean.  Furthermore,  it  is  reasonable  to  conjecture  that 
the  marginal  PDF  of  each  of  the  quadrature  components  is  symmetric 
about  the  mean  if  the  phase  is  independent  of  the  envelope  (or 
RCS)  and  the  phase  distribution  is  symmetric  about  its  mean  (but 
not  necessarily  uniform) .  In  summary,  it  is  important  to  keep  in 
mind  the  potential  impact  of  any  assumption  made  with  respect  to 
the  distribution  of  the  phase  and/or  the  dependence  (or  lack 
thereof)  between  the  phase  and  envelope  processes. 


25 


For  the  temporal  distributions  in  Table  2-2  there  are  several 
open  questions  of  relevance  to  the  approach  pursued  in  this  Phase 
I.  Two  of  those  questions  are: 

(1)  Is  the  phase  distribution  always  symmetric  about 
the  mean,  irrespective  of  the  corresponding 
envelope  PDF? 

(2)  Is  the  marginal  PDF  of  the  quadrature  components 
symmetric  about  the  mean  irrespective  of  the 
envelope  PDF? 

These  questions  are  difficult  to  answer.  One  reason  is  the  fact 
that  the  temporal  PDF  models  in  Table  2-2  are  models  of 
representation,  except  for  the  Rice  envelope  model.  Nevertheless, 
it  is  relevant  to  pursue  the  answers  to  these  and  other  related 
questions  because  the  unqualified  feasibility  of  using  third-order 
cumulants  in  a  model-based  detection  algorithm  is  determined  by 
such  issues. 

The  discussion  thus  far  has  centered  on  the  temporal  PDF 
models  in  Table  2-2  due  to  the  fact  that  cumulant-based  techniques 
are  better  suited  to  such  cases.  It  is  appropriate  to  consider 
the  cases  involving  spatial  fluctuations  of  the  clutter 
parameters.  Non-Gaussian  spatial  fluctuations  arise  as  a  result 
of  non-homogeneities  of  the  ground  (land  or  water)  clutter  over  an 
area  which  covers  a  large  number  of  radar  resolution  (range- 
azimuth)  cells.  Clutter  temporal  fluctuations  in  each  resolution 
cell  are  Gaussian-distributed,  but  the  average  RCS  fluctuates  from 
cell  to  cell  in  a  random  manner.  This  condition  is  modeled  using 
a  compound  distribution  to  describe  the  overall  clutter  return 
covering  a  large  number  of  resolution  cells,  in  contrast  to  a 
local  neighborhood  of  resolution  cells  (Lewinski,  1983;  Nathanson, 
1991) .  A  compound  distribution  model  for  the  clutter  envelope  is 
a  PDF  model  of  the  form 
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where  P(a)  is  the  overall  spatial  clutter  envelope  PDF,  P(al  r)  is  the 
clutter  envelope  PDF  at  each  individual  resolution  cell,  and  P(r)  is 
the  PDF  of  the  cell-to-cell  mean  envelope  fluctuations.  The  sea 
clutter  envelope  model  listed  in  the  second  row  of  Table  2-2  is  an 
example  of  a  spatial  variation  model  generated  from  a  compound 
distribution.  The  overall  clutter  envelope  is  K-dis tributed,  the 
conditional  PDF  is  Rayleigh,  and  the  cell-to-cell  mean  envelope 
fluctuations  behave  according  to  the  (un-normalized)  chi 
distribution  (Watts,  1987) . 

The  K-distributed  PDF  model  discussed  above  is  a  physical 
model.  Such  a  scenario  can  be  modeled  effectively  using  SIRPs, 
for  which  a  wealth  of  results  is  available  (see,  for  example, 
Rangaswamy  [1992]  and  the  references  therein) .  In  particular,  for 
a  K-distributed  SIRP  the  marginal  PDF  of  each  of  the  quadrature 
components  is  the  so-called  generalized  Laplace  PDF,  which  is 
symmetric  about  the  mean.  The  spatial  multivariate  K-dis tribution 
PDF  model  is  well  defined  for  complex-valued  as  well  as  real¬ 
valued  formulations . 

In  summary,  key  aspects  of  the  probabilistic  description  (the 
marginal  PDF  of  the  quadrature  components,  in  particular)  of  the 
non-Gaussian  clutter  types  listed  in  Table  2-2  are  unavailable, 
and  may  remain  so  for  a  long  time.  Nevertheless,  it  can  be  stated 
that  detection  techniques  based  on  higher-order  cumulants  can  be 
applied  to  the  multichannel  detection  problem  for  the  following 
two  reasons.  First,  even-ordered  cumulants  are  non-zero  for  all 
non-Gaussian  PDFs,  whether  symmetric  or  asymmetric.  This  implies 


that  the  methodology  developed  in  Phase  I  can  be  applied  using 
fourth-order  cumulants .  Second,  the  SSC  detection  methodology  can 
be  applied  in  the  cases  where  the  quadrature  components  are 
symmetrically-distributed  by  introducing  a  pre-processing  step,  as 
discussed  in  Section  3.3.  In  essence,  symmetrically-distributed 
clutter  is  transformed  by  a  zero-memory  nonlinearity  to  induce 
asymmetry  in  the  PDF  of  the  quadrature  components,  while 
introducing  tolerable  distortion  to  the  temporal  correlation 
function  (or,  equivalently,  the  frequency-domain  spectrum) .  The 
pre-processing  step  considered  in  Section  3.3  is  the  logarithmic 
function,  which  is  applied  to  the  clutter  process  prior  to 
processing.  The  degree  of  distortion  depends  on  the  dynamic  range 
and  the  frequency  content  of  the  quadrature  components,  as  well  as 
on  the  operating  region  in  the  logarithm  function.  After  the 
transformation  the  clutter  process  is  characterized  by  a  different 
univariate  PDF  with  a  larger-valued  coefficient  of  skewness  (see 
Section  5.1). 

2.1.3  Non-Gaussian  Target  Probabilistic  Models 

In  a  typical  surveillance  scenario  the  target's  physical 
dimensions  are  small  in  relation  to  the  radar  footprint;  thus,  in 
general  the  target  is  illuminated  completely  by  the  radar  beam  and 
is  present  in  the  return  from  only  one  resolution  cell.  This  is 
in  contrast  with  the  clutter,  which  has  a  large  spatial  extent 
that  covers  many  resolution  cells.  As  a  result,  the  target  RCS 
fluctuations  in  a  surveillance  scenario  of  the  type  considered 
herein  arise  only  in  the  temporal  sense. 

During  Phase  I  SSC  searched  the  literature  for  non-Gaussian 
target  models.  This  activity  led  to  the  realization  that  the 
radar  return  from  many  surveillance  targets  is  characterized  by 
non-Gaussian  statistics.  As  in  the  case  of  clutter,  the  most 
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common  probabilistic  models  for  targets  are  in  the  form  of  a 
univariate  PDF  which  characterizes  the  instantaneous  RCS 
fluctuations,  and  a  temporal  correlation  function.  Table  2-3 
lists  several  univariate  PDF  models  for  the  temporal  RCS 
fluctuations  and  the  target  types  which  can  be  represented  by  each 
PDF  model .  This  table  has  been  adapted  from  the  text  by  Nathanson 
(1991)  ,  and  enhanced  by  SSC  with  the  entry  for  the  beta 
distribution.  Notice  that  aircraft,  missiles,  and  similar  targets 
are  listed,  which  indicates  that  most  of  the  airborne  targets  of 
interest  exhibit  non-Gaussian  statistical  features. 


TEMPORAL  RCS 
FLUCTUATION  PDF 

TARGET  TYPE  DESCRIPTION 

Non-central  gamma 
(Rice  envelope) 

Targets  composed  of  one  large  reflector  and  many  small  reflectors 
(the  sum  of  the  RCS  of  all  the  small  reflectors  can  be  different  than  the 
RCS  of  the  large  reflector). 

Beta 

Complex  targets  such  as  aircraft  and  missiles;  postulated  as  a  physical 
model  as  well  as  a  representation  model  (Maffett,  1989). 

Log-normal 

Targets  with  large  mean-to-median  ratios  (>  1.44);  this  includes 
battleships,  missiles,  satellites,  and  aircraft  (at  near  broadside). 

Chi-squared  of  degree  4 
(Swerling  cases  3  and  4) 

Targets  composed  of  one  large  reflector  and  many  small  reflectors 
(the  sum  of  the  RCS  of  all  the  small  reflectors  is  equal  to  the  RCS  of 
the  large  reflector).  In  general,  the  Rice  PDF  is  a  better  model. 

Chi-squared  of  degree  2m 
(for  m  >  2) 

Short-term  (or  reduced-aspect  viewing)  statistics  of  several  types  of 
aircraft.  The  degree  is  an  inverse  function  of  the  level  of  the 
fluctuations  about  the  mean  RCS.  Equivalent  to  gamma  distribution 
for  non-integer  values  of  m. 

Weinstock 

This  model  is  included  in  the  chi-squared  PDF  of  degree  2m  for  0.3  < 
m  <  2.0. 

Table  2-3 .  Non-Gaussian  target  temporal  RCS  fluctuation  PDF 

models . 

Other  non-Gaussian  target  models  have  been  developed  for 
parameters  different  than  the  RCS.  For  example,  Novak,  Sechtin, 
and  Cardullo  (1989)  have  demonstrated  that  the  polarimetric  return 
fluctuations  from  complex  ground  targets  such  as  tanks  and  trucks 
can  be  modeled  accurately  and  effectively  using  the  K- 
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distribution.  Their  approach  is  analogous  to  the  modeling  of 
spatial  variation  using  the  compound  distribution  model  as 
summarized  in  Section  2.1.2,  with  target  aspect  angle  variability 
playing  the  role  of  spatial  ocean  surface  variability. 

The  issues  and  open  questions  identified  in  Section  2.1.2  are 
just  as  relevant  to  radar  return  from  targets.  In  particular,  the 
PDF  models  for  targets  are  deficient  in  the  same  manner  as  for 
clutter: 

multivariate  quadrature  component  PDF  for  each  of  the 
non-Gaussian  models  (except  Rice  envelope) ; 
phase  model  information  (joint  PDF  with  envelope; 
marginal  phase  distribution;  symmetry  of  marginal  phase 
PDF)  ; 

univariate  marginal  PDF  of  the  quadrature  components . 

In  the  same  manner,  candidate  resolution  approaches  identified  for 
clutter  are  relevant  also  for  targets.  This  includes  the  use  of  a 
zero-memory  nonlinear  transformation  to  skew  the  marginal  PDF  of 
the  quadrature  components . 

2.1.4  Surveillance  Radar  Scenarios  With  Non-Gaussian  Signals 

Based  on  the  clutter  and  target  PDF  models  presented  in 
Sections  2.1.2  and  2.1.3,  several  military  and  non-military 
surveillance  scenarios  and/or  detection  problems  for  radar  arrays 
can  be  identified  wherein  at  least  one  non-Gaussian  component  is 
part  of  the  radar  return.  A  variety  of  generic  scenarios  and 
detection  problems  wherein  non-Gaussian  components  arise  are 
presented  in  Table  2-4  for  military  applications,  and  in  Table  2-5 
for  non-military  applications.  The  first  two  scenarios  in  Table 
2-4  cover  the  type  of  detection  problems  that  arise  in  programs 
such  as  the  Space-Time  Processing  Program  at  RL . 
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RADAR  ARRAY 

PLATFORM 

SURVEILLANCE 

COVERAGE 

TARGETS  OF  INTEREST 

High-altitude  aircraft 

Ocean 

•  Missiles  (cruise) 

•  Aircraft  (bombers,  fighters) 

High-altitude  aircraft 

Ground 

•  Aircraft  (bombers,  fighters) 

•  Missiles  (cruise) 

Large  ships 

Airborne  threats 

•  Missiles  (cruise,  short-range) 

•  Aircraft  (bombers,  fighters) 

Aircraft 

Ground 

•  Land  vehicles  (tanks,  trucks,  etc.) 

•  Aircraft  (fighters,  helicopters) 

Table  2-4.  Candidate  military  surveillance  radar  problems 
involving  non-Gaussian  clutter  and/or  targets. 


RADAR  ARRAY 

PLATFORM 

FUNCTION 

TARGETS  OF  INTEREST 

Ground-based 

system 

Air  traffic  control 

Aircraft,  including  helicopters 

Ground-based 

system 

Weather  radar 

Weather  phenomenology 

Table  2-5.  Candidate  non-military  surveillance  radar  problems 
involving  non-Gaussian  clutter  and/or  targets. 

The  radar  return  (x(n)}  for  each  of  the  radar  surveillance 
scenarios  summarized  in  Tables  2-4  and  2-5  can  be  represented  by 
at  least  one  of  the  first  three  cases  listed  next  (presented  here 
assuming  the  alternative  hypothesis  is  true) : 

(2-20a)  x(n)  :  NGS  +  NGC  +  GWN 
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(2-20b) 

x(n): 

NGS 

+ 

GC 

+ 

GWN 

( 2  —  2  0c ) 

x(n): 

GS 

+ 

NGC 

+ 

GWN 

(2-20d) 

x(n)  : 

GS 

+ 

GC 

+ 

GWN 

where  GS  and  NGS  denote  Gaussian  and  non-Gaussian  target  signal, 
respectively;  GC  and  NGC  denote  Gaussian  and  non-Gaussian 
clutter,  respectively;  and  GWN  denotes  Gaussian  white  noise.  The 
methodologies  discussed  in  this  report  address  Case  (2-20a) .  It 
appears  that  Cases  (2-20b)  and  (2-20c)  may  be  handled  via 
modifications  and/or  extensions  to  the  approach  defined  herein. 
This  will  be  investigated  in  Phase  II.  HOS-based  identification 
and  detection  methods  are  inappropriate  for  Case  (2-20d) ,  because 
all  three  components  are  Gaussian-dis tributed .  Use  of  the  log 
transformation  may  allow  the  handling  of  Case  (2-20d)  also. 

2 . 2  Time  Series  Models 


The  time  series  model  class  is  adopted  in  this  program  to 
represent  the  target  and  clutter  components  in  the  radar  return, 
{x(n)}.  Specifically,  it  is  presummed  herein  that  Cases  (2-20a)-(2- 
20c)  can  be  expressed  in  the  form 

(2-21)  x(n)  =  y(n)  +  w(n) 

where  the  process  {y(n)}  is  the  output  of  a  time  series  system  that 
represents  the  target  and/or  clutter  components,  depending  upon 
which  hypothesis  is  valid.  For  the  data  in  one  CPI, 


2-22a) 

H0: 

y(n) 

=  C(n) 

0  <  n  <  N-1 

2~22b) 

Hi: 

y(n) 

=  s(n)  +  c(n) 

0  <  n  <  N-1 
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MA,  AR,  and  ARMA  models  are  adopted  herein  to  generate  the  process 
{y(n)},  with  emphasis  on  MA  and  AR  models  because  their  use  leads  to 
simpler  algorithms  with  reduced  analytical  and  computational 
complexity.  Haykin  and  co-workers  have  shown  that  these  models 
provide  accurate  representation  of  single-channel  radar  signals, 
specially  clutter  return  (Haykin,  Currie,  and  Kesler  [1982]; 
Metford  and  Haykin  [1985]).  Michels  (1991)  has  demonstrated  the 
value  of  model-based  multichannel  detection  using  AR  models  and 
second-order  statistics  for  Gaussian  conditions,  Case  (2-20d) . 

Moving-Average  Process.  A  J-dimensional ,  zero-mean,  stationary, 
complex-valued  sequence  {y(n)}  is  referred  to  as  a  moving- average 
vector  process  of  order  M  ,  and  is  denoted  as  MA(M),  if  it 
satisfies  a  vector  recursion  of  the  form 

M 

(2-23)  y(n)  =  X  Bku(n-k) 

k=0 


where  M  is  a  non-negative  integer,  and  {Bk  I  k  =  0,  1 . M}  are  JxJ, 

time-invariant,  complex-valued  matrices.  Equation  (2-23)  is  a 
linear  system,  with  {u(n)}  as  the  input,  {y(n)}  as  the  output,  and  {Bk} 

as  the  system  parameters.  In  the  context  of  this  program  the 
output  sequence  is  non-Gaussian;  however,  an  MA(M)  process  can  be 
Gaussian-distributed. 

Auto-Regressive  Process.  A  J-dimensional,  zero-mean,  stationary, 
complex-valued  sequence  {y(n)}  is  referred  to  as  an  auto- regressive 
vector  process  of  order  L.  and  is  denoted ‘as  AR(L),  if  it  satisfies 
a  vector  recursion  of  the  form 

L 

(2-24)  y(n)  =  -  X  Aky(n-k)  +  u(n) 

k=1 
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where  L  is  a  non-negative  integer,  and  {Ak  I  k  =  1 , 2,  .  .  .  ,  L}  are  JxJ, 

time-invariant,  complex-valued  matrices.  Equation  (2-24)  is  a 
linear  system,  with  {u(n)}  as  the  input,  {y(n)}  as  the  output,  and  {Ak} 

as  the  system  parameters.  The  output  sequence  is  non-Gaussian . 

Auto-Regressive  Moving-Average  Process.  A  J-dimensional ,  zero- 
mean,  stationary,  complex-valued  sequence  {y(n)}  is  referred  to  as 
an  auto-regressive  moving-average  vector  process  of  order  (L.M), 
and  is  denoted  as  ARMA(L,M),  if  it  satisfies  a  vector  recursion  of 
the  form 

L  M 

(2-25)  y(n)  =  -  £  A^(n-k)  +  £  B«u(n-k) 

k=1  k=0 


here  L  and  M  are  non-negative  integers  with  L>M,  and  {Aklk=1,2,  .. 
.  ,  L}  and  {Bk  I  k  =  0,  1 ,  .  .  .  ,  M}  are  JxJ,  time-invariant,  complex-valued 
matrices.  Equation  (2-25)  is  a  linear  system,  with  {u(n)}  as  the 
input ,  {y(n)}  as  the  output,  and  {Ak},  {Bk}  as  the  system  parameters. 
The  output  sequence  is  non-Gaussian. 

Characteristics  and  Properties  of  Time  Series  Models .  The 
material  presented  in  the  remainder  of  this  subsection  applies  to 
all  three  time  series  models;  distinctions  specific  to  individual 
models  are  noted.  The  input  process  {u(n)}  is  a  J-dimensional, 
zero-mean,  stationary,  complex-valued,  non-Gaussian,  white  noise 
sequence  with  finite  cumulants  at  least  up  to  sixth-order  (this 
last  assumption  is  required  in  Section  4.1)  .  The  second-order 
statistics  of  the  input  noise  are  established  as  (for  all  n 
because  the  sequence  is  stationary) , 

(2-26a)  Ruu(0)  =  =  E[u(n)u»] 
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(2-26b)  Ruu(m)  =  E[u(n)uH(n-m)]  =  [0] 


m  *  0 


u1(n)uH(n-m) 
u  u0(n)uH(n-m) 

(2-26c)  Ruu(m)  =  E[u(n)®uH(n-m)]  =  E  2  =  [0]  m  ^  0 

_uJ(n)uH(n-m)_ 

here  Zu  is  a  JxJ  Hermitian  matrix.  The  symbol  0  denotes  the 
Kronecker  matrix  product  (see,  for  example,  Pease-  [1965]),  which 
is  defined  implicitly  in  Equation  (2-26c).  Kronecker  product 
notation  simplifies  analytical  expressions  in  many  cases,  such  as 
those  involving  third-order  cumulants .  Since  the  mean  of  the 
input  sequence  is  zero,  Equations  (2-26)  define  the  covariance 
sequence  of  the  input  process,  and  matrix  £u  is  the  covariance  of 

u(n) .  The  matrix  elements  of  the  covariance  sequence  are  the 
second-order  cumulants  of  the  process  (Mendel,  1991).  Equations 
(2-26)  can  be  combined  into  one  equation  of  the  form, 

(2-27)  RUu(m)  =  E[u(n)0uH(n-m)]  =  £u  8(m)  Vm 

where  8(m)  is  the  one-dimensional  ( 1 — D )  discrete  impulse  function, 

'1  m  =  0 

(2-28)  8(m)  =  • 

0  m  *  0 

The  third-order  statistics  of  the  input  noise  sequence  are  given 
as  (for  all  n  since  {u(n)}  is  stationary) 

k  =  1, 2, . . .  ,  J 

ID-,,  m2  *  0 

k  =  1, 2 . J 


(2-29a)  Cu;k(0,0)  =  Tk  =  E[u(n)uH(n)uk*(n)] 

(2-29b)  Cu;k(m1,m2)  =  E[u(n)uH(n-m1)uk*(n-m2)]  =  [0] 


where  {rk  I  k  =  i,  2 . J}  are  JxJ ,  time- invariant ,  complex-valued 

matrices.  Using  the  two-dimensional  (2-D)  discrete  impulse 
function  (which  attains  unity  value  only  at  the  origin  of  the  2-D 

plane) ,  these  two  sets  of  equations  can  be  combined  into  one  set 

of  equations  of  the  form 

(2-30)  Cu;k(m-|,m2)  =  E[u(n)uH(n-mi)uk*(n-m2)]  =  rk5(m1,m2)  Vm1,m2 

k=  1, 2 . J 


Finally,  the  set  of  Equations  (2-30)  can  be  compacted  into  a 
single  equation  using  Kronecker  product  notation;  this  results  in 

(2-3la)  Cu(m1(m2)  =  E[u*(n-m2)®u(n)®uH(n-m1)]  Vmi,m2 


(2-31b)  Cu(mlPm2) 


SOn^  m2) 


V  nri!,  m2 


Since  the  input  sequence  is  zero-mean,  Equations  (2-31)  constitute 
the  third-order  cumulants  of  the  input  process,  and  matrices  {rk} 

are  the  non-zero  third-order  cumulants  of  the  input  process 
(Mendel ,  1991 ) . 

In  general,  the  complex-valued  cumulants  do  not  have  a 
special  structure.  However,  when  the  input  noise  is  spatially- 
independent  the  covariance  and  the  third-order  cumulants  {rk}  become 

diagonal  matrices.  The  covariance  of  a  spatially-independent 
noise  vector  u(n)  is 
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0  ...  0 


[2-32: 


Iu  =  diagLcfJ  = 


of 


0 


of 


0  0 


0 


j  J 


with  of  representing  the  variance  of  the  ith  element  of  u(n).  The 
kth  cumulant  matrix  of  a  spatially- independent  noise  vector  u(n)  is 
diagonal  also,  with  only  one  non-zero  element  located  at  the  kth 
diagonal  position: 


(2-33a) 


kth  column 

i 


0  ■ 

■  0 

0 

0  ■ 

•  •  0 

0  • 

■  0 

0 

0  ■ 

•  0 

0  • 

•  0 

\ 

0  • 

.  0 

0  • 

■  0 

0 

0  • 

•  0 

_0 

•  0 

0 

0  ■ 

•  0  _ 

<-  kth  row 


(2-33b)  E[Ui(n)Uj*(n)uk*(n)]  =  Yk  5jjk 


here  the  complex-valued  scalar  Yk  denotes  the  third-order  cumulant 
of  uk(n),  and  5yk  denotes  the  Kronecker  delta,  which  is  defined  as 


(2-34) 


f  1 


5(m)  - 


0 


i  =  j  =  k 
otherwise 


Thus,  the  third-order  statistics  of  a  J-dimensional ,  zero-mean, 

complex-valued,  spatially-independent ,  non-Gaussian ,  white  noise 
sequence  are  the  complex-valued  scalars  {Yk  1  k  =  1 , 2,  .  .  .  ,  J} . 
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For  the  models  in  Equations  (2-23 )- (2-25 )  the  input  noise  and 
the  output  sequence  are  both  non-Gaussian.  However,  it  is 
important  to  note  that  in  general  the  PDF  of  the  output  belongs  to 
a  different  family  than  the  PDF  family  of  the  input.  In  all  three 
models  the  input  process  is  transformed  into  the  output  process, 
and  these  transformations  alter  the  PDF.  Furthermore,  the  output 
process  tends  to  become  Gaussian  for  rather  general  conditions 
dictated  by  the  central  limit  theorem  (Thomas  [1969]  discusses 
one  of  the  various  versions  of  the  theorem)  .  Consider,  for 
example,  the  MA  model  as  in  Equation  (2-23) .  For  an  MA(M)  system, 
the  output  process  tends  to  be  Gaussian-distributed  for  large 
values  of  the  system  model  order,  M,  and/or  large  values  of  the 
dimension  of  the  input  noise  vector,  J.  If  a  sufficiently  large 
number  of  elements  of  the  parameter  matrices  have  non-negligible 
magnitude,  then  Gaussianity  may  be  approximated  at  smaller  values 
for  M  and/or  J.  Analogous  comments  apply  to  the  AR  and  ARMA 
models . 

The  fact  that  (for  the  three  time  series  models)  the  linear 
operations  on  the  input  sequence  transform  the  PDF  of  the  input 
sequence  creates  significant  difficulties  in  the  generation  of 
simulated  data  for  algorithm  performance  evaluation  (Section  5.1). 
This  fact  also  has  significant  implications  in  the  formulation  of 
real-time  multichannel  detection  architectures. 

An  important  construct  associated  with  the  time  series 
systems  (2-23 )- (2-25 )  is  the  system  transfer  function.  For  a  time 
series  the  system  transfer  function  is  determined  using  the  Z- 
transform  (Oppenheim  and  Schafer,  1975).  The  discrete- time  system 
transfer  functions  for  the  systems  (2-23 )- (2-25)  are  given  as, 

(2-35)  Tma(z)  =  B(z) 
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(2-36) 


Tar(z)  =  A'^z) 
Tarma(z)  =  A-^z)  B(z) 


(2-37) 

where  A(z)  and  B(z)  are  the  following  matrix  polynomials  in  Z, 

L 

(2-38a)  A(z)  =  XA|^Z'k 

k=0 

(2  —  3  8b )  A0  =  lj 

M 

(2-39)  B(z)  =  X  Bkz'k 

k=0 


Definition  (2-38b)  follows  trivially  from  Equations  (2-24) -  (2-25)  . 
The  determinants  of  the  two  polynomial  matrices,  A(z)  and  B(z),  are 
defined  as 

(2-40)  a(z)  =  IA(z)l 

(2-41)  (3(z)  =  IB(z)l 

The  properties  of  these  scalar  polynomials  determine,  in  part,  the 
behaviour  of  the  time  series .  The  matrix  pair  {A(z),  B(z)}  is 

referred  to  as  a  matrix  polynomial  representation  for  the  system 
(including  the  MA  and  AR  cases,  where  A(z)  =  lj  and  B(z)  =  lj , 

respectively) .  The  pair  (A(z),  B(z)}  is  referred  to  also  as  a  matrix 
fraction  description  (MFD)  for  the  system. 

The  conditions  satisfied  by  the  time  series  models  are  listed 
in  Table  2-6.  With  respect  to  Table  2-6,  Assumptions  (a),  (b) , 

and  (d)  insure  that  each  of  the  scalar  polynomials  a(z)  and  (3(z)  is 
of  maximum  possible  order:  a(z)  is  of  order  JL,  and  (3(z)  is  of  order 
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JM.  A  time  series  system  (MA,  AR,  or  ARMA)  for  which  Assumptions 
(c)  and/or  (e)  are  satisfied  is  referred  to  as  a  minimum  phase 
system,  and  is  said  to  have  a  minimum  phase  transfer  function  (all 
the  finite  multivariable  system  poles  and  zeros  of  a  minimum  phase 
system  are  inside  the  unit  circle) .  A  minimum  phase  system  (MA, 
AR,  or  ARMA)  and  its  inverse  system  are  both  asymptotically  stable 
(Oppenheim  and  Schafer,  1975) .  In  model-based  detection  using 
innovations  the  identified  model  must  be  minimum  phase  in  order 
for  the  whitening  filter  (model  inverse)  to  be  well  behaved  (the 
output  of  an  unstable  system  grows  unbounded) .  Assumption  (f) 
insures  that  the  ARMA  system  is  irreducible ;  that  is,  there  are  no 
pole-zero  cancelations  in  the  transfer  function  TARMA(z) 

(Rosenbrock,  1970) .  A  system  which  admits  pole-zero  cancelations 
can  be  represented  as  an  irreducible  system  of  lower  model  order. 


ASSUMPTIONS  FOR  THE  TIME  SERIES  MODELS 

MR  and  RRMR 

(a) 

B0  Has  Full  Rank 

(b) 

Bm  Has  Full  Rank 

(c) 

Zeros  of  P(z)  fire  Inside  the  Unit  Circle 

RR  and  RRMR 

(d) 

AL  Has  Full  Rank 

(e) 

Zeros  of  a(z)  Rre  Inside  the  Unit  Circle 

RRMR 

(f) 

{A(z),  B(z)}  is  relatiuely  prime  from  the  left 

Table  2-6.  Assumptions  for  the  time  series  models. 
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2 . 3  Third-Order  Cumulants 


The  cumulants  of  the  output  of  a  time  series  model  possess 
several  general  properties  of  interest,  and  satisfy  constraints 
due  to  the  structure  imposed  by  the  system  model.  Those 
constraints  are  the  basis  for  cumulant-based  model  parameter 
identification  algorithms. 

2.3.1  Definitions,  Properties,  and  Other  Issues 

The  third-order  cumulants  of  a  J-dimens ional ,  zero-mean, 
stationary,  complex-valued,  discrete- time  process  {x(n)}  are  the 
following  complex-valued  matrices : 

(2-42a)  Cx;k(m1,m2)  =  E[x(n)xH(n-m-|)xk*(n-m2)]  Vm1,m2 

k  =  1, 2 . J 

(2-42b)  C^m^nrig)  =  E[x*(n-m2)®x(n)®xH(n-m-|)]  V  rrii,  m2 

here  Xk(n)  is  the  kth  element  of  x(n) ,  matrix  Cx.k(mi,m2)  is  JxJ ,  and 
matrix  Cx(m-|,m2)  is  J2xJ.  The  second  expression  is  the  compact  form 
of  the  definition  using  Kronecker  matrix  product  notation  (the  J 
cumulants  for  a  fixed  lag  pair  are  concatenated  into  a  J2xJ 
matrix) .  This  definition  has  been  introduced  already  in  Section 
2.2  for  two  special  cases  (a  white  noise  vector  sequence,  and  a 
spatially- independent  white  noise  vector  sequence) .  The  cumulants 
{Cx(m1,m2)}  can  be  interpreted  as  a  2-D  matrix  sequence,  analogous 
to  the  1-D  covariance  matrix  seouence .  A  1-D  slice  of  {Cx(m1,m2)} 
is  a  subset  of  the  third-order  cumulants  defined  by  a  set  of  lag 
pairs  constrained  to  satisfy  a  fixed  linear  relation  between  the 
lag  indices.  For  example,  one  possible  "vertical  slice"  is 
defined  by  m-j  =  0  and  m2  =  {  .  .  .  -1 , 0,  1 ,  .  .  .  };  a  possible  "diagonal  slice" 
is  defined  by  m-|  =  m2  and  m2  =  {  .  .  .  -1 , 0,  1 ,  .  .  .  } . 
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The  (i,j)th  element  of  the  kth  cumulant  matrix  Cx;k(mi,m2)  is  a 
complex-valued  scalar  denoted  as 

(2-43)  cx;ijk(m1,m2)  =  E[xj(n)xj*(n-m1)xk*(n-m2)] 

Lower-case  notation  is  used  for  the  cumulant  in  this  equation  to 
emphasize  that  it  is  a  scalar.  In  contrast  to  second-order 
cumulants  (covariances)  for  complex-valued  vector  processes,  the 
third-order  cumulants  defined  in  Equation  (2-42)  do  not  satisfy 
symmetry  conditions.  However,  third-order  cumulants  for  scalar 
and/or  real-valued  processes  do  satisfy  symmetry  conditions,  as 
discussed  thoroughly  in  the  literature  (real-valued  scalar: 
Raghuveer  and  Nikias  [1985];  complex-valued  scalar:  Jouny  and 
Moses  [1992];  real-valued  vector:  Raghuveer  [1986]).  For  example, 
based  on  definition  (2-43)  it  is  easy  to  verify  that  for  a 
complex-valued  scalar  sequence  {h(n)}, 

(2-44)  =  ch(m1,m2) 

Other  equivalences  can  be  established  for  complex-valued  scalar 
sequences .  A  greater  number  of  equivalences  can  be  established 
for  real-valued  sequences  (vector  and  scalar) . 

Alternative  definitions  for  third-order  cumulants  are 
possible  in  the  complex-valued  vector  case,  with  the  number  and 
placement  of  conjugate  operators  as  the  major  factors  which 
determine  the  definitions.  Each  valid  definition  can  lead  to  a 
different  set  of  cumulant  constraints  as  a  result  of  the  structure 
imposed  by  the  time  series  models .  SSC  is  unaware  of  open 
literature  work  involving  cumulants  for  complex-valued  vector 
processes,  and  the  approach  and  results  presented  herein  address 
only  a  subset  of  all  possible  issues  that  can  be  identified.  The 
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complex-valued  scalar  case,  however,  has  been  considered  by  Jouny 
and  Moses  (1992)  and  by  Jelonnek  and  Kammeyer  (1992) . 

Jouny  and  Moses  (1992)  have  shown  that  the  third-order 
cumulants  must  be  defined  appropriately  in  the  context  of  AR 
models  for  a  process  defined  as  the  sum  of  harmonically-related 
complex  sinusoids  with  uniformly-distributed  phase  and  quadratic 
phase  coupling  (three  frequencies  are  harmonically-related  when 
the  highest  frequency  is  the  sum  of  the  other  two;  three 
harmonically-related  frequencies  have  quadratic  phase  coupling 
when  the  phase  associated  with  the  highest  frequency  is  equal  to 
the  sum  of  the  phases  associated  with  each  of  the  other  two 
frequencies) .  As  indicated  by  Raghuveer  and  Nikias  (1985),  phase- 
coupled  harmonically-related  complex  sinusoids  have  non-zero 
third-order  cumulants,  whereas  the  third-order  cumulants  are  zero 
if  the  phase  coupling  is  zero.  Jouny  and  Moses  (1992)  demonstrate 
that  for  such  processes  some  of  the  alternative  definitions  for 
third-order  cumulants  lead  to  incorrect  AR  parameter  estimates,  in 
the  sense  that  the  scalar  AR  model  transfer  function  does  not 
generate  the  true  bispectrum,  which  is  known  to  be  the  sum  of 
pairs  of  impulse  functions  (the  bispectrum  of  a  scalar  discrete¬ 
time  process  is  the  2-D  Fourier  transform  of  the  third-order 
cumulants) . 

The  Jouny-Moses  third-order  cumulants  for  a  complex-valued 
scalar  process  {h(n)}  are  defined  as 

(2-45)  ch(mi,m2)  =  E[{h*(n)h(n-m1)h(n-m2)}  +  {h(n)h*(n-m1)h(n-m2)} 

+  {h(n)h(n-m-|)h*(n-m2)}]  Vm1,m2 

A  different  font  element  is  used  in  this  definition  to  emphasize 
the  fact  that  Ch(m1tm2)  is  not  equal  to  Ch(mi,m2).  Notice  that 
cumulant  Ch(m1-m2)  is  defined  as  the  sum  of  three  terms,  and  that 
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each  of  the  individual  terms  can  be  used  by  itself  to  establish 
another  alternative  definition  for  complex-valued  scalar  third- 
order  cumulants.  The  extension  of  the  Jouny-Moses  definition  to 
the  vector  case  is  obvious.  However,  the  vector  version  of 
Equation  (2-45)  leads  to  intractable  expressions  in  the  derivation 
of  recursions  for  the  cumulants  of  the  output  of  a  time  series 
model  as  a  function  of  the  cumulants  of  the  input.  Furthermore, 
satisf actory  model  parameter  estimation  results  have  been  obtained 
to  date  for  non-sinusoidal  time  series  processes  using  Definition 
(2-42)  .  This  indicates  that  the  problems  observed  by  Jouny  and 
Moses  (1992)  may  be  restricted  to  sums  of  pairs  of  coupled 
complex-valued  sinusoids  with  uniformly-distributed  phase  (the 
bandwidth  of  a  complex  sinusoid  is  infinitesimally  narrow) . 
Therefore,  the  cumulant  definition  in  Equation  (2-42)  is  adopted 
for  this  program  (a  deviation  is  indulged  upon  in  Section  2.3.2 
for  the  sake  of  notational  simplicity  in  a  special  context) . 

Cumulants  posses  several  properties  which  motivate  the 
approach  of  this  Phase  I  program.  Mendel  (1991)  presents  an 
extensive  list  of  properties  of  cumulants,  with  their  associated 
proofs .  He  also  points  out  that  cumulants  can  be  treated  as  an 
operator,  much  like  the  expectation  operator.  The  properties  most 
relevant  to  the  work  herein  are  listed  in  Table  2-7. 

It  is  appropriate  to  consider  the  implications  of  these 
properties  to  radar  surveillance  problems.  Properties  (A) -(C)  in 
Table  2-7  provide  the  major  justification  for  the  application  of 
higher-order  cumulants  to  the  surveillance  radar  problem  in 
scenarios  with  non-Gaussian  clutter  and/or  targets .  Based  on 
these  properties  and  on  the  formulation  established  in  the 
preceding  sections,  it  follows  that  the  third-order  cumulants  of 
the  process  {x(n)}  are  equal  to  the  cumulants  of  the  process  which 
represents  the  non-Gaussian  clutter  and/or  target, 
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(2-46) 


Cx(m1,m2)  =  CyCm^mg)  +  CJm^rri;,)  =  Cy(m1,m2) 


This  equation  states  that  the  Gaussian-distributed  noise  is 
eliminated  from  formulations  and  computations  based  on  third-order 
cumulants .  Elimination  of  the  Gaussian-distributed  noise  effects 
is  only  approximated  in  practice  because  the  presence  of  Gaussian 
noise  generates  inaccuracies  in  the  estimates  of  the  cumulants, 
and  in  the  subsequent  estimates  of  model  parameters.  Furthermore, 
even  if  the  true  cumulants  are  utilized  to  generate  the  filter 
parameters,  noise  residuals  are  present  in  the  filter  output. 


SELECTED  PROPERTIES  OF  HIGHER-ORDER  CUMULANTS 

(R)  Cumulants  of  a  sum  of  independent  random  uariates 
equal  the  sum  of  the  cumulants  of  the  indiuidual 
uariates. 

(B)  Cumulants  of  a  set  of  random  uariates  of  which  at 
least  two  uariates  are  pairwise  independent  are 
equal  to  zero. 

(C)  Cumulants  of  a  Gaussian-distributed  scalar  random 
sequence  are  equal  to  zero. 

(D)  Odd-ordered  cumulants  of  a  scalar  random  sequence 
distributed  according  to  a  symmetric  PDF  are  equal 
to  zero. 

(E)  Cumulants  of  a  set  of  scalar  random  uariates  are 
inuariant  to  permutations  in  the  order  of  the 
arguments. 


Table  2-7.  Selected  properties  of  higher-order  cumulants. 
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Properties  (C)  and  (D)  have  been  referred  to  in  the  preceding 
sections.  These  properties  restrict  the  utilization  of  cumulants 
to  cases  involving  non-Gaussian  clutter  and/or  targets  distributed 
according  to  asymmetric  PDFs.  Lastly,  Property  (E)  is  the  basis 
for  the  equivalence  in  Equation  (2-44).  This  property  is  useful 
also  in  the  generation  of  cumulant  recursions,  as  indicated  next 
in  Section  2.3.2. 

2.3.2  Third-Order  Cumulant  Relations  for  Time  Series  Systems 

Third-order  cumulants  of  the  output  sequence  of  time  series 
systems  satisfy  algebraic  expressions  which  involve  the  input 
sequence  parameters  and  the  model  parameters  in  a  manner  analogous 
to  the  relations  satisfied  by  output  covariances  (second-order 
cumulants) .  Those  relations  are  developed  in  this  section  for  MA, 
AR,  and  ARMA  systems.  The  material  presented  in  this  section 
extends  the  multichannel  work  of  Giannakis,  Inouye,  and  Mendel 
(1989)  and  Raghuveer  (1986)  to  the  complex-valued  case,  with  the 
appropriate  modifications  to  account  for  the  specific  cumulant 
definition  adopted  in  this  program.  Equation  (2-42). 

MA  System.  Consider  the  MA(M)  process  relation,  Equation  (2-23), 
and  let  bj.k  and  by.k  denote  the  jth  column  and  the  (i,j)th  element  of 
matrix  Bk,  respectively.  That  is, 


"  b11;k 

b12;k 

b1J;k 

~  t— 1;k  — 2;k  ' 

II 

-of 

b21;k 

b22;k 

b2J;k 

bJ1;k 

bJ2;k 

bJJ;k 

It  follows  from  Equations  (2-23)  and  (2-47)  that  the  jth  element  of 
vector  y(n)  can  be  expressed  as 
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M  M  J  * 

1 2-48 )  Yj(n)  =  X  b"s^n'S)  =  X  X  bij;sUi(n'S) 

s=0  ’  s=0  i=1 


Given  these  definitions,  the  jth  third-order  cumulant  sequence  of 
{y(n)}  is  obtained  as 


'I  M  \  f  M  W  M  J 

(2-49)  Cy;j(m1,m2)  =  E  X  Bk  U(n-k)  £  uH(n-mrr)Br  £  X  bij;su*(n-m2-s) 

LI  k=0  I  \r=0  I  ls=0  i=1 


In  this  equation  the  order  of  the  expectation  and  finite 
summations  can  be  interchanged.  Carrying  out  this  interchange  and 
re-ordering  terms  results  in  the  expression 

M  M  M  J  * 

(2-50)  =  E  X  X  Z  Bk  E[u(n-k)uH(n-m1-r)ui  (n-m2-s)J  Br  bij;s 

k=0  r=0  s=0  i=1 


The  expected  value  inside  the  summations  is  recognized  as  the 
third-order  cumulants  of  the  input  noise.  Thus,  substitution  of 
Equation  (2-30)  into  Equation  (2-50)  leads  to 


M  M  M  J 


(2-5la)  Cy;j(m1,m2)  =  XZSSBf  fj  SCm^r-k.m^s-k)  Br  bij;£ 


k=0  r=0  s=0  i=1 


(2 -5  lb)  Cylm^)  =  [0] 


0  <  m-i,  m2  <  M 


m-i  >  M;  m2  >  M 


Notice  that  Equation  (2-51a)  is  defined  for  a  limited  range  of 
values  of  the  two  lags,  m-j  and  m2,  and  that  the  cumulants  are  zero 

for  all  lag  value  pairs  where  at  least  one  of  the  two  lags  is 
greater  than  the  model  order,  M.  This  is  a  result  of  the  temporal 
independence  (whiteness)  of  the  input  sequence  and  of  the  finite 
order  of  the  MA  process  (finite  summation  limit) . 
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The  remaining  steps  in  the  derivation  apply  only  to  the  non¬ 
zero  cumulants,  Equation  (2-51a) .  The  2-D  discrete  delta  function 
in  Equation  (2-51a)  assumes  the  value  of  unity  only  when  both  of 
its  arguments  are  equal  to  zero  simultaneously.  This  implies  that 
in  the  summation  over  r  all  the  terms  are  eliminated  except  the 
term  corresponding  to  r  =  k  - 17^  ;  similarly,  in  the  summation  over  S 
all  the  terms  are  eliminated  except  the  term  corresponding  to  S  =  k 
-  m2.  Thus,  Equation  (2-51)  reduces  to 

M  J 

(2-52a)  Cy;j(m-|,m2)  =  X  X  Bk  r.  Bk-mi  b.j^-rr^  0<m1lm2  <M 

k=m  i=1 


( 2 —52b)  Cy^m^iTVj) 


M 

-  s 


B 


k=m 


I  t  bii; 

i=1 


ij;k-m2 


B 


k-mi 


0  <  m2  <  M 


(2-52c)  m  =  max(m-,,m2) 

Notice  that  the  lower  limit  in  the  summation  for  k  is  now 
determined  according  to  Equation  (2-52c).  This  reflects  the  fact 
that  the  MA  model  parameter  matrices,  {Br},  are  defined  only  for 
indices  in  the  range  0  <  r  <  M,  with  r  =  k  -  rrij  and  1=1,2.  The 
expressions  in  Equations  (2-52a)  and  (2-52b)  are  equivalent; 
however,  the  second  expression  is  preferred  because  it  collects 
all  the  terms  of  the  summation  for  i.  Notice  that  each  of  the 
terms  in  the  i  summation  is  a  matrix;  it  is  convenient  to  represent 
this  sum  of  matrices  as  the  following  matrix 

J 

(2-53)  Bj(k-m2)  =  X  ri bij;k-m2  m2  <  M;  0<k<m;  1  <  j  <  J 

i=1 


The  letter  "B"  is  selected  to  represent  the  summation  in  this 
definition  in  order  to  maintain  association  with  the  MA  model 
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parameters  {bjj;k};  a  different  font  type  is  selected  to  emphasize 
that  matrices  {Bj(k-m2)}  are  distinct  from  the  MA  model  matrices  {Bk}- 

Substituting  Definition  (2-53)  into  Equation  (2-52b)  leads  to 
the  desired  final  expression, 

M 

(2-54a)  Cy.jim^mg)  =  X  Bk  Bj(k'm2)  Bk-m,  0<m1,m2<M;  1<j<J 

k=m 

(2-54b)  m  =  maxOn^rr^) 

(2-54c)  Cy;j(mi,m2)  =  [0]  m-i  >  M;  m2>M;  1  <  j  <  J 

Notice  that  Equations  (2-51b)  and  (2-52c)  are  repeated  as 
Equations  (2-54c)  and  (2-54b) ,  respectively,  to  emphasize  the  link 
that  exists  among  these  expressions. 

Equation  (2-54a)  relates  the  model  matrix  parameters  {Bk}  and 
the  matrix  cumulants  of  the  input  noise  process  {rk}  to  the  matrix 
cumulants  of  the  MA(M)  model  output  {Cy;j(mi  ,m2)}  in  a  nonlinear 
manner.  This  equation  is  similar  to  the  result  derived  by 
Giannakis,  Inouye ,  and  Mendel  (1989)  for  the  real-valued 
multichannel  case,  and  is  the  basis  for  the  MA  model  parameter 
identification  algorithm  presented  in  Section  4.2.  Equation  (2- 
54c)  provides  a  guideline  for  MA  model  order  selection  in  the 
algorithm. 

AR  Svstem:  Standard  Cumulant  Formulation.  Consider  the  AR(L) 
process  relation,  Equation  (2-24),  and  let  a.j;k  and  denote  the 
jth  column  and  the  (i,j) th  element  of  matrix  Ak,  respectively.  That 
is, 
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al1;k 

al2;k 

a1J;k 

Ak  _  C^1;k  ^ 2;k  ’  ' 

•  -J;J  “ 

a21;k 

a22;k 

a2J;k 

aJ1  ;k 

aJ2;k 

aJJ;k 

As  in  the  MA  case,  the  definition  introduced  in  Equation  (2-55) 
can  be  used  to  express  the  jth  element  of  the  AR  vector  y(n)  as 

L_  L  J 

(2-56)  y-(n)  =  -  X  a^sy(n-s)  +  uj(n)  =  -  X  X  ajjsyi(n's)  +  uj(n) 

s=1  ’  s=1  i=1 


Substitution  of  Equation  (2-24)  for  y(n-m1)  into  the  jth  third-order 
cumulant  definition,  Equation  (2-42a) ,  results  in  the  following 
expression. 


(2-57)  Cy;j(mi,m2)  =  E 


Y(n)  {-  X  yH(n-m1-k)  Ak)y*(n-m2) 

k=1 


+  E 


y(n)  u^n-m^  y*(n-m2) 


The  order  of  the  expectation  and  finite  summation  on  the  first 
term  on  the  right-hand-side  of  this  equation  can  be  interchanged. 
Also,  the  second  term  on  the  right-hand-side  can  be  expanded  by 
substitution  of  Equation  (2-24).  These  manipulations  lead  to  the 
expression 


(2-58)  Cy.j(m1,m2)  =  -  X  E  y(n)yVmrk)  y,  (n-m2)  A 


k=1 


+  E 


X  Ak  y(n-k) )  u^n-iTii)  y* (n-m2) 

k=1 


+  E 


u(n)  ^(n-m^y:  (n-m2) 
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The  expectation  inside  the  summation  on  the  first  term  on  the 
right-hand-side  of  this  equation  is  the  jth  third-order  cumulant  at 
lags  (m-|+k,m2) ,  and  the  order  of  the  expectation  and  finite 

summation  can  be  interchanged  on  the  second  term.  Also,  Equation 
(2-56)  can  be  substituted  for  yj(n-m2)  on  the  third  term.  These 

operations  result  in 

L 

(2-59)  Cy;j(mi,m2)  =  -  X  Cy;j(m1+k,m2)  Ak 

k=1 

L 

-  X  Ak  E[  y(n-k)  uH(n-m1)  y*(n-m2) 

k=1 

+  E  u(n)  uH(n-m1)  X  X  aij;sy*(n"m2-s)  j 
+  E  u(n)  uH(n-m.,)  u*(n-m2) 

The  order  of  the  expectation  and  finite  summations  can  be 
interchanged  on  the  third  term  on  the  right-hand-side  of  this 
equation.  Also,  the  last  term  is  recognized  as  the  jth  third-order 
cumulant  of  the  input  noise  sequence.  Carrying  out  these 
manipulations  leads  to  (after  re-ordering  the  terms) 

L 

(2-60)  Cy;j(m.|,m2)  =  -  X,  Cy;j(m1+k,m2)  Ak  +  r^m^nr,) 

k=1 
L 

-  X  AkE  y(n-k)  u^n-rrij)  y*(n-m2) 

k=1 

L  J  *  1 

-  X  X  E[y(n)uH(n-m1)y‘(n-m2-s)Jaij;s 

s=1  i=1 


1  <  j  <  J 


Equation  (2-60)  is  a  general  expression  which  is  valid  for  all 
values  of  the  lags  ITl-(  and  m2 .  Further  substitutions  of  the  AR 

expression  into  Equation  (2-60)  do  not  lead  to  a  simpler  form. 
However,  an  alternative  approach  results  in  a  simplified  relation, 
as  indicated  next. 

The  number  of  terms  on  the  right-hand-side  of  Equation  (2-60) 
can  be  reduced  by  selecting  a  subset  of  all  possible  lag  values 
such  that  Property  (B)  in  Table  2-7  applies.  Specifically,  for 
lags  m1  =  m2  =  0  the  last  two  terms  in  the  right-hand-side  are  equal 
to  zero,  and  for  lags  in  the  range  0  <  rr^  <  m2  the  last  three  terms 
in  the  right-hand-side  are  equal  to  zero.  That  is, 

L 

(2-61a)  Cy-j(0,0)  =  -  y,  Cy;j(k,0)  Ak  +  Tj  1  <j<J 

k=i 

L 

(2-6lb)  Cy;j(mi,m2)  =  -  ^  Cy.^rr^+k.rr^)  Ak  0  <  m1  <  m2;  1  <  j  <  J 

k=1 


These  equations  can  be  written  in  a  form  which  emphasizes  their 
recursive  structure  and  is  related  to  covariance  formulations.  In 
order  to  do  so  it  is  convenient  to  first  define  a  JLxJ  block  column 
matrix  fl  as 


(2-62) 


fl  = 


*1 

A0 


Now  it  is  possible  to  write  Equations  (2-61a)  in  the  form  of  a 
single  block  equation  as 


(2-63) 


Cy(0,0)  +  [Cy(1,0)  Cy(2,0)  •••  Cy(L,0)]H  =  Cu(0,0) 


where  the  definitions  in  Equations  (2-31)  and  (2-42b)  have  been 
used.  The  block  row  matrix  which  pre-multiplies  H  in  Equation  (2- 
63)  is  dimensioned  J  xJL.  Consider  now  Equations  (2-61b)  for  a  1-D 
slice  defined  by  m1  =  0  and  m2  =  {1 , 2,  .  .  .  ,  S},  with  S>L.  For  this  1-D 
slice  Equations  (2-61b)  can  be  expressed  in  expanded  form  as 


(2-64) 


_Cy(1,1) 

Cy(2,1)  • 

Cy(L,1)_ 

"V 

_Cy(0,1)" 

Cy(1,2) 

Cy(2,2)  • 

Cy(L,2) 

A2 

=  - 

Cy(0,2) 

_Cy(1,S) 

Cy(2,S)  ■ 

••  Cy(L,S) 

aJ 

_Cy(0,S) 

Now  define  a  J2xJL  block  row  matrix  Cy;1L;  a  J2SxJL  block  matrix 
Cy:s,L '  anc^  a  J2SxJ  block  column  matrix  Cy.gj  as 


(2-65)  Cy:1iL  =  [Cy(1 ,0)  Cy(2,0) 


Cy(L,0)] 


(2-66) 


"Cy(1,1)  Cy(2,1) 

Cy(1,2)  Cy(2,2) 


Cy(L,1) 

Cy(L,2) 


Cy(I.S)  Cy(2,S)  Cy(L,S) 


Cy(0,1) 

Cy(0,2) 

Cy(0,S) 


Combining  the  definitions  in  Equations  (2-62 )- (2-67 ) ,  it  follows 
that 
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(2-68) 


Cy:S,L  A  ~  "  Cy:Sli 


(2-69)  CU(0,0)  =  Cy(0,0)  +  Cyl  |_  fl 

Equation  (2-68)  is  referred  to  as  the  third-order  cumulant  normal 
equation ,  and  the  J2SxJL  matrix  Cy;siL  is  the  normal  matrix.  This 
normal  matrix  does  not  have  a  particular  structure. 

Equation  (2-68)  is  used  to  identify  the  AR  matrix  parameters 
(Section  4.3) .  Given  the  AR(L)  matrix  coefficients,  Equation  (2- 
68)  is  used  to  generate  an  estimate  of  the  input  noise  third-order 
cumulants .  AR  model  recursions  of  this  type  are  discussed  by 
Mendel  (1991)  and  others  for  the  real-valued  scalar  case. 

The  form  of  Equations  (2-68)  and  (2-69)  results  from  the 
approach  adopted  for  the  derivation,  and  from  the  specific  1-D 
slice  selected  to  constrain  the  lag  values.  Different  choices 
(for  the  approach  and/or  the  1-D  slice)  result  in  different 
expressions  with  distinct  analytical  and  numerical  properties.  An 
option  of  interest  leads  to  a  normal  matrix  with  block  Toeplitz 
structure,  as  summarized  next. 


Recursive  equations  with  a  structure  different  from  that  of 
Equations  (2-61)  can  be  generated  by  changing  the  steps  in  the 
derivation  presented  above.  Specifically,  consider  the  third- 
order  cumulant  definition,  Equation  (2-42a) ,  and  substitute  the  AR 
system  expression  for  y(n),  instead  of  for  y(n-m-|)  as  in  Equation  (2- 

57).  After  a  few  manipulations  similar  to  those  carried  above, 
the  following  expressions  are  obtained: 


(2-70a) 


Cy;j(0,0) 


L 

=  X 

k=l 


Ak  Cy;j(-k,-k)  +  r 


1  <j  <  J 
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(2—7  Ob )  Cy;j(m1,m2)  =  -  £  Cy;j(m1-k,m2-k) 

k=1 


m1(  m2  >  0;  1  <  j  <  J 


These  equations  are  analogous  to  Equations  (2-61) ,  but  exhibit  a 
different  structure  (it  is  important  to  note  that  in  Equation  (2- 
70b)  it  suffices  that  either  m-|  >  0  or  m2  >  0)  .  As  before,  the  J 

Equations  (2-70a)  can  be  combined  into  a  single  expression.  Also, 
for  Equations  (2-70b)  select  the  diagonal  1-D  slice  defined  by  m-|  = 
m2  =  {1,  2,  .  .  .  ,  L},  and  combine  all  J  equations  into  a  single 
expression.  Finally,  define  the  following  three  block  matrices: 


(2-71) 


(2-72) 


Cy:L,1  - 


-y:L,L  “ 


Cy(-I.-I)  " 
Cy(-2,-2) 

Cy(-L,-L)  _ 
Cv(0,0) 


rCy;l(-1,-1) 

Cy;j(-1,-1)_ 

Cy;1  (-2,-2)  ..- 

Cy;j(-2,-2) 

[Cy^-L.-L)  ... 

Cy;j(-L,-L). 

Cy(1,1) 

•  Cy(L-I.L-l) 

)  Cy(0,0) 

•  Cy(L-2,L-2) 

+1)  Cy(-L+2,-L+2) 

Cy(0,0) 

(2-73) 


Cy:i,L  =  [Cy(1,1)  Cy(2,2)  Cy(L,L)] 


Cy.|_ti  is  a  JLxJ2  block  column  matrix,  Cy:L>L  is  a  JLxJ2L  block  matrix, 

o 

and  Cy. i  |_  is  a  JxJ  L  block  row  matrix.  Given  these  definitions, 
Equations  (2-70)  are  expressed  compactly  as 


[2-1  A) 


C,Hl.l  A 


:2-75)  cu(o,o)  =  [ r,  r2 


rj]  =  cy(0.0)  +  nHcy:UI 
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where  fl  is  as  defined  in  Equation  (2-62) .  The  Hermitian  operator 
is  used  in  Equation  (2-74)  so  that  this  equation  can  be  compared 
directly  with  Equation  (2-68)  .  Notice  that  matrix  Cy:Li|_  is  block 

Toeplitz  and  block  square.  Also  notice  that  the  arrangement  of 
the  matrix  elements  as  a  block  row  in  Cy(m,m)  and  in  cu(0,0)  differs 

from  the  definition  used  heretofore  (block  column  arrangement) . 
This  notational  deviation  is  introduced  to  obtain  compact  block 
matrices  in  Equations  (2-74)  and  (2-75),  instead  of  sparse, 
higher-dimensioned  matrices .  The  fact  that  an  alternative 
notation  is  preferable  when  a  different  derivation  approach  is 
pursued  provides  another  indication  of  the  richness  of  options 
available  when  HOS  for  vector  processes  are  involved. 

Equations  (2-74)  and  (2-75)  contain  equivalent  information  as 
Equations  (2-68)  and  (2-69)  with  S  =  L.  However,  it  is  likely  that 
each  of  these  two  pairs  of  equations  will  lead  to  different 
results  in  the  practical  solution  of  modeling,  identification,  and 
related  problems.  An  extended  version  of  Equation  (2-74)  is 
obtained  by  using  S  >  L  index  pairs  (mi,m2)  in  the  diagonal  1-D 

slice.  The  resulting  matrix  corresponding  to  Equation  (2-72)  has 
block  Toeplitz  structure,  but  is  not  block  square. 

AR  System:  Raahuveer  Cumulant  Formulation.  Raghuveer  (1986)  has 
proposed  an  AR  model  recursion  based  on  a  subset  of  the  full  set 
of  third-order  cumulants  for  the  real-valued  multichannel  case. 
SSC  extended  Raghuveer ' s  formulation  to  the  complex-valued 
multichannel  case.  Raghuveer ' s  formulation  requires  fewer 
computations  than  the  normal  equations  approach  (since  fewer 
scalar  cumulants  are  estimated) ,  and  the  matrix  generated  by  the 
recursion,  although  not  a  normal  matrix,  has  Toeplitz  structure. 
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Raghuveer ' s  formulation  is  based  on  a  "reduced"  cumulant 
matrix.  Given  the  AR(L)  system  output  sequence  {y(n)},  define  a  JxJ 
matrix  Y(n)  as 


(2-76) 


Y(n)  =  diag[yi(n)]  = 


y1  (n)  o  ... 

0  y2(n)  ... 


0 

0 


L  0  0  ...  yj(n)J 


Let  Cy(m)  denote  a  JxJ  matrix  of  third-order  cumulants  defined  as 
(2-77)  Cy(m)  =  E[y(n)yH(n-m)Y*(n-m)]  V  m 

The  elements  of  this  matrix  are  a  subset  of  the  third-order 
cumulants  of  the  process  {y(n)}.  Thus,  matrix  Cy(m)  is  not  equal  to 

the  third-order  cumulants,  Cy(m1,m2),  as  defined  in  Equation  (2-42); 

3 

in  fact,  matrix  Cy(m1,m2)  has  J  elements  (third-order  scalar 

2 

moments)  whereas  matrix  Cy(m)  has  only  J  elements.  Furthermore, 
matrix  Cy(m)  is  different  than  matrix  Cy(m1,m2)  restricted  to  the  1- 
D  slice  defined  by  m1  =  0  and  m2  =  m . 


For  an  AR(L)  process  {y(n)},  cumulant  matrix  Cy(m)  restricted  to 
lag  values  in  the  set  1  <m<L  results  in 


(2-78) 


Cy(m) 


E 


k=1  / 


+  E  u(n)  yH(n-m)  Y*(n-m) 


1  <  m  <  L 


In  this  equation  the  order  of  expectation  and  finite  summation  can 
be  interchanged  on  the  first  term  to  the  right  of  the  equal  sign. 
Also,  for  the  selected  range  of  lags  the  second  term  on  the  right- 
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hand-side  is  equal  to  zero  (past  outputs  are  independent  of  future 
inputs) .  Thus,  Equations  (2-78)  reduce  to 


(2-79) 


1  <  m  <  L 


Since  the  AR  process  is  stationary,  each  of  the  three  factors 
inside  the  expectation  operator  can  be  shifted  forward  in  time  by  k 
integer  steps;  the  result  is  recognized  as  Cy(m-k) .  Equations  (2- 

79)  become 

L 

(2-80)  Cy(m)  =  -  X  Ak  Cy(m-k)  1<m<L 

k=1 


In  expanded  form,  Equations  (2-80)'  are  expressed  as 


Cy(0) 

Cy(1)  ■ 

••  Cy(L-l) 

Cy(-1) 

Cy(0)  • 

Cy(L-2) 

Cy(1)  Cy(2)  ■  •  ■  Cy(L)  ]  =  -  flH 

'  Cy(-L+1) 

Cy(-L+2) 

•  Cy(0) 

where  R  is  the  JLxJ  block  matrix  defined  in  Equation  (2-62)  .  Now 
define  a  JLxJL  block  matrix  Cypp  and  a  JxJL  block  row  matrix  Cy.-jp 
as  (bold  type  is  used  to  emphasize  that  these  matrices  are 
distinct  from  their  counterparts  in  Equations  (2-68)  and  (2-69), 


Cy(0) 

Cy(1) 

••  Cy(L-l)" 

Cy(-1) 

Cy(0)  • 

■■  Cy(L-2) 

Cy(-L+1) 

Cy(-L+2)  .. 

Cy(0)  _ 

(2-83)  Cy;  1  L  -  Cy(1)  Cy(2)  Cy(L) 


58 


By  applying  the  definitions  in  Equations  (2-62)  and  (2-82 )- (2-83 ) , 
the  following  compact  expression  is  obtained  for  Equation  (2-81) , 


(2-84)  c"LiLn  =  -CyH:1,L 


The  Hermitian  operator  has  been  applied  in  Equation  (2-84)  so  that 
this  expression  can  be  compared  directly  with  Equations  (2-68)  and 
(2-74)  .  The  JLxJL  block  matrix  Cy.^L  has  block  Toeplitz  structure. 


ARMA  System:  Standard  Cumulant  Formulation.  Consider  the  output 
process ,  {y(n)},  of  an  ARMA(L,M)  system.  Equation  (2-25),  and  the 
third-order  cumulant  definition  in  Equation  (2-42).  Substitution 
of  the  ARMA  definition  for  y(n)  into  Equation  (2-42a)  results  in 


(2-85) 


Cy;j(m1,m2)  =  E 


-i<y(n-k)j 

k=1  ) 


+  E 


Bg  u(n-s) 


yH(n-m1)  y*(n-m2) 
Jy^n-m^y^n-mg) 


The  order  of  the  expectation  and  finite  summation  operations  can 
be  interchanged  in  each  of  the  two  terms  on  the  right-hand-side  of 
this  equation.  This  leads  to 


[2-86)  Cy;j(m1,m2)  =  -  £  A^E[y(n-k)yH(n-m1)y*(n-m2) 

k=1 


M 


+  x  Bs  E  ^(n_s)  yH(n-mi)  V*.  (n'm2) 


s=0 


Since  both  the  input  and  output  processes  are  stationary,  the 
factors  inside  the  first  expectation  operator  can  be  shifted 
forward  in  time  by  k  integer  steps,  and  the  factors  inside  the 
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second  expectation  operator  can  be  shifted  forward  in  time  by  S 
integer  steps.  Thus,  Equations  (2-86)  become 

L  r 

(2-87)  CyjjOn^mg)  =  -  2  Ak  E  y(n)  yH(n-m1+k)  y*(n-m2+k) 

k=1 

M  H  r 

+  X  BS  ELu(n)  yH(n-m1+s)  y*(n-m2+s) 

s=0 

The  expected  value  factors  in  the  first  summation  are  recognized 
as  Cy.j(m1-k,m2-k) .  Similarly,  it  is  recognized  that  the  expected 

value  factors  in  the  second  summation  are  equal  to  zero  for  m1  >  M 
and  m2  >  0  (or  equivalently,  for  m1  >  0  and  m2  >  M)  .  Therefore,  a 
subset  of  Equations  (2-87)  is  obtained  as 

V1  H 

(2-88)  Cy;j(m1)m2)  =  -  2^  Ak  Cy;j(m1-k,m2-k)  m-,  >  M;  m2  >  0;  1<j<J 

k=i 

The  structure  of  this  equation  is  identical  to  that  of  Equation 
(2-70b) .  However,  the  range  of  applicable  values  for  m1  is  less 

in  Equation  (2-88)  than  in  Equation  (2-70b) .  Proceeding  as  in  the 
discussion  following  Equation  (2-70b) ,  the  final  result  of 
Equation  (2-75)  is  valid  for  this  case  also,  with  the  provision 
for  the  different  range  of  values  for  m1 .  Namely, 

(2-89)  CJ!LiLfl  =  -  Cy:1L  m1  >  M;  m2  >  0 

Block  matrices  Cy:L  L  and  Cy:1L  have  structure  analogous  to  the 

structure  defined  in  Equations  (2-72)  and  (2-73).  As  an  example, 

for  a  1-D  slice  defined  by  m1  =  m2  =  {L+1 ,  L+2 . 2L},  matrices  Cy.L  L 

and  Cy;1L  are  given  as 
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(2-90) 


c 


y;  l.l 


Cy(L,L)  Cy(L+1,L+1) 

Cy(L-1,L-1)  Cy(L,L) 

Cy(1,1)  Cy(2,2) 


Cy(2L-1,2L-1) 

Cy(2L-2,2L-2) 

Cy(L,L) 


(2-91) 


Cy-  1  l_  =  Cy(L+1,L+1)  Cy(L+2,L  +  2) 


Cy(2L,2L)] 


For  the  selected  1-D  slice,  the  conditions  on  m-j  and  m2  expressed 
in  Equation  (2-89)  are  satisfied  because  L>M.  Notice  that  matrix 
Cy;i,L  is  block  Toeplitz.  The  real-valued  scalar  ARMA  case  has  been 
discussed  by  Mendel  (1991)  and  others. 


ARMA  System:  Raqhuveer  Cumulant  Formulation.  Consider  the  third- 
order  cumulant  definition  in  Equation  (2-77),  applied  to  an 
ARMA(L,M)  process  {y(n)};  Equation  (2-25).  Substitution  of  the  ARMA 
definition  into  Equation  (2-77)  restricted  to  lag  values  in  the 
range  m  >  1  results  in 


(2-92)  Cv(m)  =  E 


'  X  Ak  y(n'k)  yH(n-m)  Y*(n-m) 


L  k=1 


+  E 


M 


X  Bs  u(n-s)  yH(n-m)  Y*(n-m) 


L  s=0 


m  >  1 


Following  the  approach  established  earlier,  it  is  possible  to  show 
that  Equations  (2-92)  reduce  to  the  following 

L  M  r  , 

(2-93)  Cy(m)  =  -  X  A^  Cy(m-k)  +  £  Bs  E[u(n)yH(n-m+s)  Y*(n-m+s)J  m>1 

k=1  s=0 
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Notice  that  the  second  summation  is  equal  to  zero  for  lag  values 
in  the  range  m>M  (since  past  outputs  are  independent  of  future 
inputs).  Therefore,  an  important  subset  of  Equations  (2-93)  is 

L  H 

(2-94)  Cy(m)  =  -  Ak  Cy(m-k)  m>M 

k=1 


This  equation  is  identical  to  Equation  (2-80),  except  that  the 
range  of  applicable  lag  values  for  m  in  Equation  (2-80)  is  m>1, 
instead  of  m  >  M.  Proceeding  as  in  the  discussion  following 
Equation  (2-80),  the  final  result  of  Equation  (2-84)  is  valid  for 
this  case  also,  with  the  different  range  of  values  for  m.  Namely, 

(2-95)  Cy^B  =  -  Cy.-j  L  IT1>M 

Block  matrices  Cy:(_  |_  and  Cy;1i_  have  structure  analogous  to  the 

structure  defined  in  Equations  (2-82)  and  (2-83).  For  the  case 

where  the  range  of  lag  values  is  m  =  {L+1,L+2 . 2L},  matrices  Cy:L|_ 

and  Cy;i  |_  are  given  as 


"  Cy(L) 
Cy(L-l) 

Cy(L+1) 

Cy(L) 

•••  Cy(2L-1)" 
...  Cy(2L-2) 

(2-96) 

Cy:L,L  “ 

.  Cy(1) 

Cy(2) 

•  •  •  Cy(L) 

(2-97) 

Cy:1,L  =  1 

;  Cy(L+1) 

Cy(L+2) 

...  Cy(2L)] 

Notice  that  the 

constraint  on 

m  stated  in  Equation  (3-47)  is 

satisfied 

because  L  >  M . 

Also , 

matrix  Cy;|_[_  is  block  Toeplitz. 

These  results  extend  Raghuveer ' s  formulation  to  complex-valued 
multichannel  ARMA  systems. 
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3 . 0  DETECTION  METHODOLOGY 


The  detection  case  considered  in  this  report  is  Case  (2-20a) , 
wherein  the  baseband  sequence  {x(n)}  is  composed  of  a  spatially- 
unresolved  non-Gaussian  target  embedded  in  additive  non-Gaussian 
clutter  and  additive  Gaussian  noise.  Several  issues  associated 
with  Case  (2-20a)  are  addressed  in  this  section,  and  possible 
detection  options  are  presented. 

Define  the  set  of  vectors  {x(n-1),  x(n-2), .  . . ,}  as  the  past  of  the 
sequence  {x(n)},  and  denote  as  X  (n)  the  vector  space  covering  all 
possible  linear  combinations  of  the  elements  of  the  past  of  {x(n)} . 
The  minimum  variance  estimate  (MVE)  of  x(n)  is  the  conditional 
expectation  E[x(n)IJT(n)] .  Determination  of  the  MVE  and  of  the  LLR 

detector  requires  availability  of  the  multivariate  PDF,  and  in  the 
case  of  a  non-Gaussian  sequence  the  MVE  and  the  LLR  detector  may 
be  nonlinear  functions  of  the  data.  As  mentioned  in  Sections 
2.1.2  and  2.1.3,  the  multivariate  PDF  is  unavailable  in  most  non- 
Gaussian  cases;  also,  nonlinear  functions  can  be  difficult  to 
implement.  Thus,  suboptimal  linear  approximations  to  the  MVE  and 
to  the  LLR  detector  were  adopted  in  Phase  I. 

Let  x(nin-1  ;Hj)  for  i  =  0,  1  ,  denote  a  linear  estimate  of  x(n)  based 
on  X'(n)  and  conditioned  on  hypothesis  Hj  being  true.  The  linear 
estimates  referred  to  are  as  defined  in  Section  3.5  using  model 
parameters  estimated  as  described  in  Section  4.0.  Now  define  a 
conditional  pseudo-innovations  sequence,  {e(nlHj)}  for  i  =  0,  1  ,  as 

(3-1)  e(nlHj)  =  x(n)  -  x(nln-1  ;Hj)  i  =  0,  1 

The  conditional  pseudo-innovations  is  a  zero-mean  sequence  because 
both  the  baseband  process  and  its  linear  estimate  are  zero-mean. 
This  sequence  is  not  a  true  innovations  sequence  because  the 
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estimate  x(nln-1;Hj)  as  determined  in  Section  3.5  is  not  the  MVE,  and 
because  the  noise  sequence  {w(n)}  is  non-zero  (if  w(n)  =  0  for  all  n, 
then  the  estimate  generated  as  in  Section  3.5  is  a  linear  MVE). 
For  these  same  reasons  a  detector  which  uses  the  pseudo¬ 
innovations  (3-1)  is  not  the  LLR  detector.  In  place  of  the  true, 
but  unavailable,  LLR  detector  for  the  non-Gaussian  case,  the 
approximate,  suboptimal  detector  proposed  by  Metford  (1984)  can  be 
used.  The  pseudo-innovations  of  Equation  (3-1)  is  used  in  such  a 
configuration,  as  discussed  in  Section  3.4.  In  a  minor  abuse  of 
notation,  e(nlH0)  and  e(nIH-,)  are  used  herein  to  denote  the  output  of 

the  null  hypothesis  filter  and  the  alternative  hypothesis  filter, 
respectively.  At  each  resolution  cell  only  the  output  of  one  of 
the  two  filters  is  a  valid  pseudo-innovations. 

3 . 1  Detector  Architecture  Configurations 

Two  basic  configurations  can  be  defined  for  detection 
architectures:  (a)  adaptive  on-line,  and  (b)  off-line.  The 

adaptive  on-line  configuration  is  illustrated  in  Figure  3-1.  In 
this  configuration  the  parameters  of  the  two  filters  (one  filter 
corresponding  to  each  hypothesis)  and  covariance  matrices  required 
for  detection  rule  calculations  are  determined  in  real-time  by 
processing  the  multichannel  baseband  sequence  for  at  least  one 
CPI.  Filter  order  may  be  estimated  on-line  also  to  provide 
complete  adaptability.  Order  determination  is,  however,  a 
difficult  task.  Alternatively,  the  orders  of  the  two  filters  may 
be  pre-stored  in  the  radar  processor  memory.  The  order  of  the 
alternative  hypothesis  filter  is  always  higher  than  the  order  of 
the  null  hypothesis  filter  because  the  radar  return  contains  more 
information  (the  target  component)  when  the  alternative  hypothesis 
is  true.  An  on-line  configuration  provides  the  most  adaptability, 
but  presents  a  large  computational  burden  even  if  filter  order  is 
pre-stored  as  a  function  of  scenario  conditions. 
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Higher-order  cumulants  can  be  used  to  adapt  other  aspects  of 
the  radar  processor.  As  an  example,  estimates  of  higher-order 
cumulants  can  be  used  to  assess  the  extent  of  deviation  of  the 
radar  return  from  a  Gaussian-distributed  sequence.  Such  knowledge 
can  be  used  to  select  the  detection  rule  most  appropriate  for  the 
identified  conditions. 

Pseudo- 


Figure  3-1.  Multichannel  detector  configuration  with  on-line 
adaptive  filter  parameter  identification. 

The  second  configuration  is  illustrated  in  Figure  3-2.  The 
major  differences  between  the  two  configurations  involve  the 
central  block  in  each  of  the  two  diagrams.  In  this  configuration 
each  filter  is  designed  off-line  using  simulated  and/or  true  radar 
data,  and  the  resulting  filter  designs  are  implemented  in  the 
radar  processor.  Multiple  filter  pairs  may  be  designed  off-line, 
with  each  filter  pair  tuned  to  a  distinct  set  of  operational 
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conditions.  This  provides  a  partial  degree  of  adaptability.  The 
multichannel  baseband  sequence  is  used  to  estimate  covariance 
matrices  required  for  detection  rule  calculations,  and  to  select 
the  filter  pair  which  best  matches  the  scenario  conditions  if  more 
than  one  filter  pair  is  pre-stored  in  the  radar  processor  memory. 
Careful  off-line  design  of  the  filters  can  lead  to  acceptable 
performance . 


Pseudo- 


Figure  3-2.  Multichannel  detector  off-line  architecture 
configuration  (off-line  filter  design) . 

3 . 2  Detection  Domain 


As  discussed  in  Section  2.1.2,  a  non-Gaussian  temporal 
clutter  PDF  model  represents  the  conditions  wherein  the 
instantaneous  fluctuations  of  the  clutter  signal  at  each 
individual  resolution  cell  behave  according  to  a  non-Gaussian 
distribution.  In  turn,  a  non-Gaussian  spatial  clutter  PDF  model 
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represents  conditions  wherein  the  clutter  is  non-homogeneous  over 
a  large  physical  area  which  includes  many  resolution  cells,  but 
the  instantaneous  temporal  fluctuations  in  each  individual 
resolution  cell  behave  according  to  a  Gaussian  distribution.  A 
cumulant-based  detection  architecture  can  be  designed  for  each  of 
these  two  conditions. 

Consider  first  a  scenario  wherein  the  non-Gaussian  temporal 
clutter  PDF  model  is  valid,  and  the  filter  parameters  are 
estimated  on-line.  In  such  a  scenario,  at  each  individual 
resolution  cell  the  third-order  cumulants  of  the  multichannel 
baseband  sequence  are  non-zero  and  can  be  estimated  for  each  CPI 
using  time  averages  (recall  from  Table  2-1  that  {x(n)}  is  ergodic)  . 
In  turn,  the  estimated  cumulants  are  used  in  the  identification 
algorithms  to  estimate  the  model  parameters  of  a  whitening  filter 
tuned  to  the  particular  resolution  cell.  The  procedure  is 
repeated  for  each  resolution  cell,  and  it  is  likely  that  the 
estimated  filter  varies  from  cell  to  cell.  The  domain  in  which 
this  procedure  is  carried  out  is  a  local  domain;  more 
specifically,  it  is  a  single-cell  domain.  A  detection 
architecture  defined  for  a  single  resolution  cell  is  referred  to 
herein  as  a  single-cell  detection  architecture. 

In  an  off-line  architecture  configuration  the  whitening 
filter  design  is  done  off-line.  In  this  case  the  non-Gaussian 
characteristic  of  the  radar  return  in  individual  cells  is  taken 
into  account  in  the  estimation  of  third-order  cumulants  and  in  the 
design  of  the  filters  for  use  at  each  individual  resolution  cell. 

Consider  now  a  scenario  in  which  the  non-Gaussian  spatial 
clutter  PDF  model  is  applicable.  In  this  case  all  higher-order 
cumulants  determined  at  individual  cells  over  a  CPI  are  zero. 
However,  since  the  clutter  process  is  non-Gaussian  for  the  total 
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region,  the  ensemble  of  radar  returns  for  all  resolution  cells 
over  a  CPI  represents  a  non-Gaussian  process.  An  ensemble 
estimate  (as  oppossed  to  a  time-average  estimate)  of  third-order 
cumulants  can  be  obtained,  and  the  estimated  cumulants  can  be  used 
to  identify  the  parameters  of  a  whitening  filter  which  is  tuned  to 
all  resolution  cells  in  the  region.  This  procedure  has  a  global 
domain .  and  a  detection  architecture  defined  for  a  whole  region 
with  a  large  number  of  resolution  cells  is  referred  to  herein  as  a 
global  detection  architecture . 

A  global  detection  architecture  can  be  configured  either  as 
an  adaptive  on-line  or  an  off-line  architecture.  In  a  global 
detection  architecture  a  single  filter  pair  is  used  for  a  large 
number  of  resolution  cells.  This  is  a  simple  architecture  with 
significantly  reduced  computational  requirements  in  relation  to 
the  requirements  of  a  single-cell  detection  architecure,  specially 
for  an  off-line  configuration.  The  detection  architecture 
proposed  by  Rangaswamy,  Weiner,  and  Michels  (1993)  for  SIRPs  is  a 
global  detection  architecture.  It  is  different,  however,  from  the 
cumulant -based  global  detection  architecture  introduced  herein. 

3 . 3  Pre-Processing  Option 

As  stated  in  Sections  2.1.2  and  2.1.3,  it  is  unknown  whether 
the  marginal  PDF  of  the  quadrature  components  is  symmetric  or 
asymmetric  about  the  mean  because  the  marginal  PDFs  are 
unavailable  for  the  non-Gaussian  temporal  clutter  and  target  PDF 
models  in  Tables  2-2  and  2-3.  Since  the  degree  of  asymmetry  of 
the  marginal  PDFs  is  an  open  question,  it  is  appropriate  to  define 
candidate  options  to  handle  each  of  the  two  possibilities.  If  the 
marginal  PDF  is  asymmetric,  then  the  third-order  cumulants  are 
non- zero  and  the  cumulant-based  procedures  defined  in  this  report 
can  be  applied  directly.  On  the  other  hand,  if  the  marginal  PDF 
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is  symmetric,  then  the  third-order  cumulants  are  zero.  This  case 
can  be  handled  via  two  different  approaches,  as  discussed  next. 

One  approach  is  to  formulate  the  procedures  in  this  report 
using  fourth-order  (or  higher-order)  cumulants  instead  of  third- 
order  cumulants.  This  approach  is  feasible  because  many  non- 
Gaussian  symmetric  PDFs  often  have  fourth-order  cumulants  with 
large  magnitude.  However,  in  the  multichannel  case  estimation  of 
fourth-order  cumulants  requires  significantly  more  computations 
than  estimation  of  third-order  cumulants.  Furthermore,  parameter 
identification  algorithms  based  on  fourth-order  cumulants  are 
computationally  more  intensive  than  equivalent  algorithms  based  on 
third-order  cumulants. 

The  second  approach  is  motivated  by  the  fact  that  the 
marginal  PDF  of  the  quadrature  components  can  be  altered  (skewed) 
significantly,  while  the  temporal  correlation  sequence  is  affected 
to  a  tolerable  degree,  by  the  application  of  some  zero-memory 
nonlinear  transformations  to  each  of  the  two  real-valued 
quadrature  components  prior  to  any  cumulant-based  processing.  One 
such  nonlinear  transformation  is  the  logarithm  function,  as 
indicated  by  Zheng,  McLaughlin,  and  Mulgrew  (1993) .  The  logarithm 
of  {x(n)}  is  computed  as  a  pre-processing  step  preceding  the  block 
diagrams  in  Figures  3-1  and  3-2.  That  is,  the  following 
operations  are  carried  out  on  the  multichannel  baseband  sequence, 


(3-2a) 

z,(n)  =  logb[x,(n)  +  a,] 

( 3 -2b) 

ZqOi)  =  logb[  XQ(n)  +  cJq] 

( 3 -2c ) 

|i  =  E[z,(n)  +  j  ZQ(n)] 

( 3 -2d) 

z(n)  =  z,(n)  +  j  ZQ(n)  -  \i 

69 


where  b  denotes  the  logarithm  base,  fi.  is  the  complex-valued  mean 
of  the  sequence  z,(n)  +  j  z^n) ,  and  a(  and  3q  are  real-valued  vectors 

determined  as 

(3-3a)  a,  =  lmin{X|(n)}l  +  c] 

( 3 -3b)  aQ  =  lmin{xQ(n)}l  +  d 

(3-3c)  C  >  1 

with  C  a  real-valued  scalar,  and  j  is  a  J-dimensional  vector  with 
all  elements  equal  to  unity.  The  function  of  the  vectors  a|  and  aQ 

is  to  bias  the  sequence  {x(n)}  so  that  the  argument  of  the  logarithm 
function  in  Equations  (3 -2a)  and  (3 -2b)  is  greater  than  or  equal 
to  unity  for  all  n  (recall  that  the  sequence  {x(n)}  has  zero  mean)  . 
In  a  hardware  implementation  the  sample  mean  of  Z|(n)+jZQ(n)  is  used 

in  place  of  the  unavailable  true  mean.  Notice  that  sequence  {z(n)} 
has  zero  mean,  as  required  for  cumulant  estimation.  If  this  pre¬ 
processing  option  is  used,  the  transformed  sequence  {z(n)}  replaces 
{x(n)}  in  the  detection  procedure.  However,  in  order  to  avoid 
confusion  (x(n)}  will  continue  to  denote  the  multichannel  baseband 
sequence  in  the  remainder  of  this  report. 

The  use  of  the  logarithmic  transformation  has  several 
features  of  interest.  The  log  transformation  is  instantaneous 
(zero-memory),  and  can  be  implemented  as  a  table  look-up  in  a 
digital  processor.  In  some  radar  systems  the  logarithm  is  an 
inherent  function  in  the  analog  portion  of  the  receiver,  applied 
to  the  radar  return  signal  in  order  to  reduce  its  dynamic  range 
(see,  for  example,  Nathanson  [1991]).  The  effect  of  the 
logarithmic  transformation  on  simulated  data  is  examined  in 
Section  5.0.  The  results  discussed  therein  indicate  that  the  log 
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transformation  can  be  used  effectively  for  data  generation  as  well 
as  a  pre-processor  in  the  cases  where  the  PDF  of  the  quadrature 
components  is  symmetric  about  the  mean. 

3 . 4  Detection  Rule 


Metford  (1984)  (see  also  Metford,  Haykin,  and  Taylor  [1982]) 
has  demonstrated  that  the  LLR  detection  rule  for  the  Gaussian 
case.  Case  (2-20d),  can  be  used  in  an  approximate  sense  as  the 
detection  rule  for  the  case  of  detecting  a  non-Gaussian  sequence 
embedded  in  additive  Gaussian  noise.  Furthermore,  acceptable 
results  can  be  obtained  for  the  cases  where  the  true  innovations 
sequence  is  unavailable  (or  equivalently,  where  only  a  suboptimal 
estimate  of  the  non-Gaussian  sequence  is  available) .  These 
conditions  are  similar  to  the  conditions  for  Case  (2-20a) ,  as 
described  in  Sections  2. 0-4.0  of  this  report.  Additionally,  SSC 
has  noticed  in  simulation-based  analyses  that  the  pseudo¬ 
innovations  generated  using  the  whitening  filters  for  time  series 
models  tend  to  be  Gaussian-distributed  as  the  model  order  and/or 
the  number  of  channels  increases.  This  is  a  manifestation  of  the 
central  limit  theorem:  a  linear  combination  of  a  large  number  of 
non-Gaussian  random  variates  tends  to  be  Gaussian-distributed.  It 
is  likely  that  the  pseudo- innovations  tend  to  be  Gaussian- 
distributed  also  when  there  is  some  model  mismatch  (the  true 
system  differs  from  the  model  adopted  for  the  filter  design) .  A 
degree  of  model  mismatch  is  present  when  real  radar  data  is 
concerned  because  the  time  series  models  are  inherently  models  of 
representation.  For  these  reasons  the  Gaussian  LLR  is  adopted  in 
this  program  as  an  approximate,  suboptimal  detection  rule  for  non- 
Gaussian  sequences. 

The  innovations -based  multichannel  LLR  detection  rule  for  the 
complex-valued  Gaussian  case  has  been  derived  by  Michels  (1991). 


Also,  Rangaswamy,  Weiner,  and  Michels  (1993)  have  demonstrated 
recently  that  the  methodology  can  be  extended  to  include  non- 
Gaussian  SIRPs  (recall  that  the  SIRP  model  is  appropriate  for  the 
spatial  non-Gaussian  clutter  conditions  as  described  in  Section 
2.1.2).  Based  on  the  above  discussion,  Michels'  multichannel  LLR 
detection  methodology  is  adopted  in  this  program,  with  the 
conditional  pseudo-innovations  used  in  place  of  the  conditional 
innovations.  For  brevity,  only  the  final  LLR  expression  is 
presented  here. 

Consider  the  conditional  pseudo- innovations  in  Equations  (3- 
1 )  ,  and  let  jQ(Hj)  denote  the  covariance  matrix  of  the  pseudo¬ 
innovations  , 

(3-4)  n(Hj)  =  E[  £(nlHj)  £H(nlHj)]  i  =  0,  1 


This  is  a  true  covariance  matrix  because  the  pseudo-innovations 
sequence  is  zero-mean.  Now  let  0(Ho,H-|)  denote  the  multichannel 

likelihood  ratio  for  Case  (2-20d) ,  the  Gaussian  case.  Then,  the 
complex-valued  multichannel  LLR  as  defined  by  Michels  (1991)  is 


(3-5) 


,  N-1 

)  =  I 

In 

Q(«o) 

1 

n=0 

L  0(Hi)|J 

-  £H(nlH, )  n'RHU  a(n!H, )] 

As  indicated  in  Figures  3-1  and  3-2,  the  LLR  is  compared  to  a 
threshold,  T,  which  is  calculated  adaptively  to  maintain  a 

constant  false  alarm  rate  (CFAR) , 

(3-6)  ln[0(Ho,H1)]  = 


select  H-, 
select  H0 
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A  candidate  CFAR  approach  with  demonstrated  good  performance 
consists  of  calculating  the  median  of  a  set  of  the  LLR  values  from 
a  number  of  adjacent  range  cells  (at  the  same  azimuth)  on  both 
sides  of  the  cell  in  question,  and  then  scaling  the  calculated 
median  value  by  a  pre-determined  constant  to  provide  the  desired 
false  alarm  rate  (Metford  and  Haykin,  1985)  . 

Alternative  expressions  for  the  LLR  can  be  generated  based  on 
factorization  of  the  pseudo-innovations  covariance  matrix  and 
spatial  whitening  of  the  pseudo-innovations  sequence.  Applicable 
covariance  factorization  methods  include  Cholesky  factorization, 
LDU  decomposition,  and  singular  value  decomposition  (SVD) ;  all 
three  lead  to  simplified  LLR  expressions.  The  first  two 
techniques  have  been  described  by  Michels  (1991)  ,  and  the  SVD 
technique  is  described  by  Roman  and  Davis  (1993) . 

3 . 5  Digital  Realizations  For  Whitening  Filters 

The  cumulant-based  model  parameter  identification  algorithms 
presented  in  Section  4 . 0  generate  the  parameters  of  a  time  series 
model  for  the  multichannel  baseband  sequence,  {x(n)} .  Given  the 
matrix  parameters  for  a  time  series  model,  the  whitening  filter 
which  generates  the  pseudo-innovations  is  obtained  as  the  inverse 
system  of  the  time  series  model.  Specifically,  the  whitening 
filter  for  a  MA(M)  system  is  an  AR(M)  f  ilter;  the  whitening  filter 
for  an  AR(L)  system  is  a  MA(L)  f  ilter;  and  the  whitening  filter  for 
an  ARMA(L,M)  system  is  an  ARMA(L,L)  filter . 

A  time  series  whitening  filter  can  be  implemented  with 
various  distinct  digital  realizations,  such  as  tapped  delay  line 
(or  direct  form),  parallel,  and  cascade.  These  three  types  of 
digital  realizations  are  minimal  (or  canonical)  because  their 
implementation  requires  the  minimum  number  of  delay  operations . 
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The  tapped  delay  line  realization  is  of  particular  interest 
because  it  can  be  generated  by  inspection  of  the  time  series 
system  expressions.  In  general,  cascade  and  parallel  digital 
filter  realizations  offer  the  best  numerical  performance 
(minimization  of  parameter  quantization  effects),  but  require 
knowledge  of  the  system  poles  and/or  zeros  (Oppenheim  and  Schafer, 
1975)  .  The  system  poles  and  zeros  are  the  roots  of  the 
polynomials  a(z)  and  P(z),  respectively,  as  defined  in  Equations  (2- 

40)  and  (2-41) .  Time  series  systems  which  satisfy  the  assumptions 
in  Table  2-6  have  JL  poles  and/or  JM  zeros.  Computation  of  the 
polynomials  a(z)  and  P(z)  and  their  roots  imposes  a  large 

computational  load  and  is  difficult  to  implement  on-line.  It  is 
feasible,  however,  in  the  off-line  case. 

In  the  remainder  of  this  section  the  conditional  notation  is 
dropped  from  the  argument  set  of  the  pseudo- innovations , 
estimates,  and  other  variates  for  notational  simplicity.  The 
figures  and  equations  presented  are  applicable  for  both  hypotheses 
with  the  stipulation  that  the  system  model  order  is  larger  for  the 
alternative  hypothesis  (the  target  forces  additional  modeling 
degrees-of-freedom) . 

3.5.1  Moving  Average  System  Whitening  Filter 

A  tapped  delay  line  realization  for  the  whitening  filter  of  a 
MA(M)  system  is  presented  in  Figure  3-3.  In  this  figure,  as  well 
as  in  the  next  two  figures,  each  of  the  blocks  labeled  "DELAY" 
represents  a  J-dimensional  column  of  scalar  delay  operators  because 
the  input  to  each  block  is  a  J-dimensional  vector.  This  is  a 
minimal  realization,  with  a  total  of  JM  delays.  The  linear 
estimate  and  the  pseudo-innovations  for  this  case  are 
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(3-7) 


x(nln-l)  =  £  B^e(n-k) 

k=1 

M 

(3-8)  e(n)  =  x(n)  -  x(nln-l)  =  -  ^  B^e(n-k)  +  x(n) 

k=1 

From  Equation  (3-8),  the  whitening  filter  is  an  AR(M)  system  with 
X(n)  as  the  input,  e(n)  as  the  output,  and  {Bk}  as  the  AR  matrix 

parameters .  The  MA  system  zeros  are  the  poles  of  the  whitening 
filter;  the  whitening  filter  does  not  have  any  finite-value  zeros. 


Figure  3-3 .  Tapped  delay  line  realization  for  whitening  filter 

corresponding  to  a  MA(M)  system. 


3.5.2  Auto-Regressive  System  Whitening  Filter 

A  tapped  delay  line  minimal  realization  (JL  delays)  for  the 
whitening  filter  of  an  AR(L)  system  is  presented  in  Figure  3-4. 
The  linear  estimate  and  the  pseudo-innovations  for  this  case  are 

L 

(3-9)  x(nln-l)  =  A^x(n-k) 

k=i 

L  L 

(3-lOa)  e(n)  =  x(n)  -  x(nln-l)  =  x(n)  +  Yj  Ak£(n’k)  =  X  Ak*(n'k) 

k=1  k=0 

(3-10b)  A0  =  lj 

Thus,  the  whitening  filter  is  a  MA(L)  system  with  x(n)  as  the  input, 
e(n)  as  the  output,  and  {Ak}  as  the  MA  matrix  parameters.  The  AR 
system  poles  are  the  zeros  of  the  whitening  filter;  all  the 
whitening  filter  poles  assume  the  value  zero. 

3.5.3  Auto-Regressive  Moving-Average  System  Whitening  Filter 

A  tapped  delay  line  minimal  realization  for  the  whitening 
filter  of  an  ARMA(L,M)  system  is  presented  in  Figure  3-5.  This 
realization  has  a  total  number  of  JL  delays.  The  linear  estimate 
and  the  pseudo- innovations  for  this  case  are 

L  M 

(3-11)  x(nln-l)  =  -  X  A^x(n-k)  +  £  B^e(n-k) 

k=1  k=1 


L  L 

(3-12a)  e(n)  =  x(n)  -  x(nln-l)  =  -  B^e(n-k)  +  ^  A^x(n-k) 

k=1  k=0 


76 


( 3 -12b) 


Ao  -  lj 


(3-12c)  Bk  =  [0] 


k  =  M+1,  . .  . ,  L 


From  Equation  (3-12),  the  whitening  filter  is  an  ARMA(L.L)  system 
with  x(n)  as  the  input,  e(n)  as  the  output,  {Ak}  as  the  MA  parameters, 
and  {Bk}  as  the  AR  parameters.  The  ARMA  system  zeros  are  the  poles 

of  the  whitening  filter,  and  the  ARMA  system  poles  are  the  zeros 
of  the  whitening  filter.  When  L>  M,  the  extra  poles  in  the 
whitening  filter  assume  the  value  zero. 


x(n) 


i 


DELAY 


x(n-L) 


-A' 


H 


Figure  3-4.  Tapped  delay  line  realization  for  whitening  filter 
corresponding  to  an  AR(L)  system. 
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4 . 0  IDENTIFICATION  ALGORITHMS 


Identification  of  the  model  parameter  matrices,  {Ak}  and/or 
{Bk},  for  the  whitening  filters  is  carried  out  based  on  the  third- 
order  cumulant  relations  for  time  series  systems  presented  in 
Section  2.3.2.  The  MA  and  AR  identification  algorithms  presented 
herein  can  be  categorized  as  linear  algebra  algorithms  because  in 
these  methods  the  key  step  in  the  estimation  of  the  unknown  matrix 
parameters  is  the  solution  of  a  system  of  equations  which  are 
linear  in  a  set  of  parameters.  In  comparison  with  optimization 
algorithms  (which  require  the  iterative  solution  to  a  nonlinear 
system  of  equations),  linear  algebra  algorithms  are  simpler  to 
implement  in  software  and/or  hardware,  and  do  not  involve 
convergence  issues  directly.  Optimization  algorithms  do  provide 
more  accurate  results  at  an  increased  computational  cost.  Thus, 
linear  algebra  algorithms  are  better  candidates  for  on-line 
implementation,  and  optimization  algorithms  are  better  candidates 
for  off-line  analyses  and  design.  The  identification  algorithms 
summarized  herein  require  the  third-order  cumulants  of  the 
multichannel  baseband  sequence;  since  the  true  cumulants  are 
unavailable,  an  estimate  is  obtained  and  used  in  the  algorithms. 

Stability  of  the  identified  model  and  the  associated 
whitening  filter  is  an  issue  of  relevance  to  all  known 
identification  algorithms  based  on  higher-order  cumulants. 
Specifically,  none  of  the  identification  algorithms  based  on 
higher-order  cumulants  can  guarantee  that  the  estimated  parameters 
correspond  to  a  stable  system.  This  is  due,  in  part,  to  the  fact 
that  third-order  cumulants  contain  phase  information  (unlike  the 
covariance  sequence),  which  allows  the  use  of  third-order 
cumulants  to  identify  minimum-phase  as  well  as  non-minimum  phase 
models  (recall  that  if  a  system  is  minimum-phase,  then  it  is 
stable,  and  its  whitening  filter  is  stable  also) .  Analyses 
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carried  out  to  date  by  SSC  indicate  that  all  the  algorithms 
considered  in  Phase  I  can  identify  non-minimum  phase  models  even 
when  a  stable  system  is  used  to  generate  the  data.  However,  this 
occurs  only  in  a  limited  number  of  cases.  SSC  has  observed  also 
that  MA  algorithms  are  more  prone  to  generate  unstable  system 
estimates  than  AR  algorithms.  Further  analyses  are  required 
during  Phase  II  in  order  to  establish  MA,  AR,  and  ARMA  algorithm 
performance  with  respect  to  stability. 

In  the  context  of  multichannel  detection  for  radar 
surveillance,  model  stability  is  an  issue  only  for  adaptive  on¬ 
line  detector  configurations.  Two  possible  approaches  to  address 
this  potential  problem  for  on-line  configurations  have  been 
identified.  The  first  approach  is  to  average  several  estimates  of 
model  parameters  obtained  over  a  homogeneous  region.  The  second 
approach  is  to  extend  the  detection  rule  formulation  to  handle 
non-causal  whitening  filters,  and  to  assign  the  unstable  poles 
and/or  zeros  to  the  non-causal  component  of  the  filter.  Stability 
is  not  an  issue  if  causality  is  not  a  requirement.  Since  in  an 
airborne  radar  processor  the  received  data  can  be  stored  at  least 
for  one  CPI,  it  is  possible  to  consider  a  non-causal  processing 
approach . 

4 . 1  Third-Order  Cumulants  Estimation 

Each  of  the  two  detector  architectures  mentioned  in  Section 

3.2  requires  a  different  cumulant  estimator.  The  first  estimator 
is  a  time  average  estimator  for  a  single-cell  domain  detector 
formulation  and  is  defined  as 

N-l 

(4-la)  Cx^m^mg)  =  — ! —  T  x(n) ^(n-m^ x*(n-m2)  j  =  1,2,  ...,J 

N  -  m  n=m  ‘ 
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(4-lb) 


N-1 

Cx(mi,m2)  =  — p —  X  x*(n-m2) ® x(n) ® xH(n-rrii) 

N  -  171  n=m 


(4-lc)  m  =  max(mi,m2) 

In  this  estimator  the  time  average  is  taken  over  one  CPI.  If 
appropriate,  it  is  possible  to  include  data  for  more  than  one  CPI. 
Estimator  (4-1)  provides  an  unbiased  estimate  of  the  multichannel 
baseband  sequence  third-order  cumulants  since  the  temporal  process 
is  assumed  to  be  ergodic  and  zero  mean  (Table  2-1) .  The  estimate 
(4-1)  converges  in  probability  to  the  true  third-order  cumulants 
because  {x(n)}  is  assumed  to  be  exponentially  stable,  and  the  model 
input  noise  process  {u(n)}  is  assumed  to  be  stationary  with  finite 
cumulants  at  least  up  to  sixth-order.  The  condition  of  finite 
cumulants  at  least  up  to  sixth-order  guarantees  that  the  first  six 
cumulants  of  the  input  noise  process  are  absolutely  summable, 
which  is  the  actual  requirement  for  convergence  in  probability 
(see  Mendel  [1991]  and  the  references  therein) . 

The  second  estimator  is  an  ensemble  estimator  for  a  global 
domain  detector  formulation  and  is  defined  as 

(4-2a)  6^(171! ,m2)  =  1  J  xk(n) x^(n-nii) x*k(n-m2)  j  =  1,2,  ...,J 

ke  "K 

(4-2b)  CxOn^rr^)  =  -^  X  Xk(n-m2) ® xk(n) ® x^n-m^ 

ke^C 


The  subscript  k  denotes  the  radar  resolution  cell  (as  in  Section 
2.1  and  Figure  2-1),  denotes  the  set  of  selected  ensemble 
resolution  cells,  and  K  denotes  the  number  of  resolution  cells 
included  in  %.  The  number  K  must  be  large  in  order  to  insure  that 

the  global  non-Gaussian  features  are  represented  well.  If  the 
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selected  resolution  cells  are  independent,  then  a  smaller  number 
of  cells  suffices.  The  time  n  is  fixed  at  an  instant  in  the  CPI 
for  each  resolution  cell  in  the  ensemble.  The  ensemble  mean  must 
be  estimated  and,  if  non-zero,  subtracted  from  the  ensemble 
samples  prior  to  estimation  of  the  cumulants .  Estimator  (4-2)  is 
an  unbiased  estimator  of  the  true  third-order  cumulants;  also,  the 
estimate  converges  in  probability  to  the  true  third-order 
cumulants  for  the  same  reasons  as  Estimator  (4-1) . 

The  following  approximations  are  adopted  implicitly  in  all 
equations  and  discussions  where  cumulants  appear  and  the  true 
cumulants  are  unavailable, 

(4-3a)  Cx^m-jjiTis)  —  Cx;j(m1lm2)  j  =  1, 2 . J 

(4-3b)  Cx(m1,m2)  ~  C^m^ma) 

That  is,  the  appropriate  estimate  (either  (4-1)  or  (4-2))  is  used 
in  place  of  the  unknown  true  third-order  cumulants . 

4 . 2  MA  Parameter  Identification 


The  MA  parameter  identification  algorithm  presented  herein  is 
a  modification  of  the  algorithm  proposed  by  Giannakis,  Inouye,  and 
Mendel  (1989)  .  The  modifications  made  by  SSC  include  extension  of 
the  algorithm  to  the  complex-valued  case,  and  simpler  expressions 
for  the  identification  equations.  In  this  report,  this  algorithm 
is  referred  to  as  the  one-slice  MA  algorithm. 

One-Slice  MA  Algorithm.  The  one-slice  identification  algorithm 
for  MA  systems  is  based  on  two  key  decisions  involving  the  third- 
order  cumulant  relations  for  MA  systems,  Equations  (2-54) .  First, 
a  specific  1-D  slice  of  Equations  (2-54a)  is  selected  so  that  the 
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summation  on  the  right-hand-side  has  only  one  non-zero  term.  The 
selected  1-D  slice  is  defined  by  0  <  m1  <  M  and  m2  =  M .  For  the  lag 

pairs  in  this  slice  the  summation  on  the  right-hand-side  of 
Equations  (2-54a)  includes  only  one  non-zero  term.  Second,  matrix 
Bm  is  taken  to  be  equal  to  the  J-dimensional  identity  matrix. 
This  does  not  imply  a  loss  of  generality  because  matrix  BM  has 
full  rank  (Table  2-6),  and  because  the  basis  in  which  the  noise 
input  {u(n)}  is  identified  is  irrelevant  with  respect  to  detection 
performance  (recall  also  that  the  MA  model  is  a  representation 
model)  . 


Consider  the  third-order  cumulant  relations  for  MA  systems 
derived  in  Section  2.3.2,  and  let  m2  =  M  in  Equations  (2-54a). 

Also,  concatenate  Equations  (2-54a)  for  j  =  1,2,  to  obtain 

(4-4)  Cx(m1,M)  =  [  Ij  ®B^]  B(0)  BM_mi  0  <  nr^  <  M 


2 

where  Cx(m-|,M)  and  B(0)  are  J  xJ  matrices  defined  as 


(4-5) 


Cx(mi,M) 


Cx;1(m1fM) 

Cx;2(mi,M) 


Cx;j(m1(M) 


0  <  m1  <  M 


(4-6)  B(0) 


B^O) 

B2{0) 


L  Bj(0)  J 


Now  let  m1  = 0  in  Equation  (4-4),  and  use  the  fact  that  BM  has  full 
rank  to  obtain 
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(4-7) 


Jj  ®bJJ]  B(0)  =  Cx(0,M)B^ 


This  expression  is  substituted  into  Equation  (4-4)  to  eliminate 
matrix  B(0),  which  results  in 

(4-8)  C^m^M)  =  Cx(0,M)B^BM.mi  0  <  m-|  <  M 

Since  the  basis  adopted  to  represent  the  input  noise  sequence  {u(n)} 
is  arbitrary  with  respect  to  detection  performance,  matrix  B^  can 

be  assumed  to  be  an  identity  matrix.  Given  this  assumption  and 
Equation  (4-8),  the  MA  matrix  parameters  are  obtained  as 

(4 -9a)  BM  =  lj 

(4-9b)  BM„m,  =  C^(0,M)  Cx(m1,M)  0  <  m-,  <  M 

where  the  dagger  (t)  denotes  the  pseudoinverse  operator. 

To  identify  the  third-order  cumulants  of  the  input  noise  the 
MA  matrix  parameters  are  estimated  first  according  to  Equation  (4- 
9)  .  Next,  let  m1  =  M  and  m2  =  0  in  Equations  (2-54a)  to  obtain 

(4-10)  CX;j(M,0)  =  Bj(M)  B0  j  =  1,  2 . J 

with  Equation  (4-9a)  taken  into  account.  From  Definition  (2-53), 

J 

(4-11)  Bj(M)  =  Xribij:M  j  =  1.  2 . J 

i=i 

Since  bjj;M  is  the  (i,j)th  element  of  BM,  it  follows  from  Equation  (4- 
9a)  that 
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(4  12)  by;M  —  5jj 

Consequently,  Equation  (4-11)  reduces  to 

(4-13)  Bj(M)  =  Tj  j  =  1, 2,  .  .  .  ,  J 

And  Equations  (4-10)  become 

(4-14)  Cx;j(M,0)  =  Tj  B0  j  =  1 , 2, .  . . ,  J 

Matrix  BQ  is  determined  using  Equation  (4-9b)  for  m-|=M,  and  is  a 
full  rank  matrix  (Table  2-6)  .  Therefore,  the  parameters  {11}  are 

identified  according  to  the  expression 

(4-15)  Tj  =  Cx;j(M,0)  B'o  =  Cx;j(M,0) [cJ(0,M)  CX(M,M)]  1  j  =  1, 2 . J 

These  cumulants  are  not  required  for  the  detection  rule. 
Computation  of  the  pseudoinverse  in  Equation  (4-9b)  and  of  the 
inverse  in  Equation  (4-15)  can  be  carried  out  via  the  SVD. 

An  alternative  multichannel  MA  identification  algorithm  has 
been  proposed  by  Tong,  Inouye,  and  Liu  (1992)  .  Their  algorithm 
converges  in  a  finite  number  of  steps,  but  is  significantly  more 
complicated  than  the  one-slice  algorithm  presented  above.  The 
Tong- Inouye-Liu  algorithm  should  be  considered  in  Phase  II, 
including  evaluation  of  its  capability  to  generate  correct 
parameter  estimates  for  minimum-phase  systems. 

Model  Order  Determination.  Model  order  for  an  MA(M)  system  can  be 

determined  based  on  Equation  (2-54c) .  Assume  that  the  true  model 
order  is  less  than  a  given  upper  bound,  denoted  as  M j .  Then 

compute  as  many  1-D  slices  as  possible  of  the  third-order 
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cumulants  Cx(m-|,m2).  At  least  two  1-D  slices  should  be  computed. 

Two  convenient  1-D  slices  are  the  horizontal  and  vertical  slices, 
defined  as  {0  <  m1  <  MT,  m2  =  0}  and  {mi=0,  0<m2<MT},  respectively. 
In  theory,  the  cumulants  for  lag  pairs  with  either  IT^  >  M  or  m2  >  M 

are  zero,  but  in  practice  all  cumulants  are  likely  to  be  non-zero. 
However,  if  good  estimates  of  the  cumulants  Cx(m1fm2)  are  available, 

then  a  discontinuity  in  the  magnitude  of  the  elements  of  the 
estimated  cumulant  matrices  can  be  expected  at  lag  pairs  with 
either  m1  =  M+1  or  m2  =  M+1  .  The  next  step  in  the  procedure  is  to 

apply  the  principle  of  a  matrix  norm  to  each  of  the  estimated 
cumulant  matrices  (a  matrix  norm  is  useful  in  assessing  relative 
magnitude  between  square  matrices  [Faddeeva,  1959]).  The  set  of 
computed  norms  is  a  2-D  function  of  m-)  and  m2 .  This  2-D  function 

is  examined  to  determine  the  points  of  discontinuities,  which 
provide  an  estimate  of  model  order.  A  candidate  matrix  norm, 
referred  to  herein  as  matrix  norm  A,  is  defined  for  a  JxJ  matrix  C 
with  elements  {Cy},  as 


(4-16) 


IICII 

A 


max 

j 


2 

Since  the  cumulant  matrices  are  rectangular  (J  xJ)  ,  matrix  norm  A 
cannot  be  applied  directly.  However,  matrix  norm  A  can  be  applied 
to  a  JxJ  matrix  computed  as 


(4-17)  C  =  CxOm^mg)  Cx(m.|,m2) 


Other  matrix  norms  are  available  (such  as  the  Euclidean  matrix 
norm  [Faddeeva,  1959]),  but  matrix  norm  A  is  attractive  because 
squaring  operations  are  not  involved  (by  Equation  (4-17),  the 
elements  of  the  cumulants  are  squared  already) . 
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4 . 3  AR  Parameter  Identification 


The  AR  parameter  identification  algorithms  summarized  herein 
are  discussed  by  Mendel  (1991)  for  the  real-valued  scalar  case, 
and  Raghuveer  (1986)  for  the  real-valued  multichannel  case.  The 
algorithms  are  based  on  the  recursive  relations  for  cumulants 
derived  by  SSC  for  the  complex-valued  multichannel  case  and 
presented  in  Section  2.3.2.  Two  distinct,  but  similar,  algorithms 
are  summarized.  The  first  algorithm  is  based  on  the  standard 
cumulant  formulation,  and  the  second  algorithm  is  based  on  the 
Raghuveer  cumulant  formulation. 

AR  Model  Identification  Using  the  Standard  Cumulant  Formulation. 
Consider  the  third-order  cumulant  recursions  presented  as 
Equations  (2-68)  and  (2-69),  and  recall  that  S > L  +  1  .  Equation 
(2-68)  can  be  viewed  as  a  collection  of  J  linear  systems  of 
equations  (corresponding  to  the  J  columns  of  fl) ,  with  each  system 
consisting  of  JSxJL  equations  in  JL  unknowns  (the  total  number  of 
scalar  unknowns  is  J  L)  .  Thus,  Equation  (2-68)  is  an  over- 
determined  system  if  the  rank  of  the  J  SxJL  matrix  Cx;S  |_  is  equal 
to  JL.  Matrix  Cx:Sil  has  full  rank  in  most  cases  if  S>L+1  (see 
Mendel  [1991]  and  the  references  therein),  as  selected  herein. 
The  linear  system  of  Equations  (2-68)  is  solved  via  the 
pseudoinverse  operator  to  result  in 

(4-18)  fl  =  -  C^.S  L  Cx:Sj1 

Using  the  matrix  parameters  obtained  via  Equation  (4-18),  the 
input  noise  third-order  cumulants  are  estimated  using  Equation  (2- 
69)  directly,  repeated  here  as 

(4-19)  Cu(0,0)  =  Cx(0,0)  +  Cx:1  L  H 
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(recall  that  cu(0.0)  is  a  block  column  matrix  with  Tj  as  the  jth 
block  element)  .  The  pseudoinverse  in  Equation  (4-18)  can  be 
implemented  using  the  SVD  or  another  robust  technique. 

The  SVD  can  be  used  also  to  determine  model  order.  As  in  the 
case  of  the  MA  system,  assume  that  an  upper  bound  for  the  true 
model  order  is  known,  and  denote  this  upper  bound  as  Ly.  Let  {Sj  I  i  = 
1,2,...,  Ly}  denote  the  singular  values  of  matrix  Cx:Sii_,  arranged  in 

decreasing  order  of  magnitude  (as  usual) .  In  ideal  circumstances, 
only  the  first  JL  singular  values  are  non-zero,  but  in  practical 
cases  all  singular  values  are  likely  to  be  non-zero.  Given  good 
estimates  of  the  cumulants,  a  significant  discontinuity  in  the 
magnitude  of  the  singular  values  can  be  expected  after  SL.  Thus, 

the  SVD-based  model  order  estimate  is  the  number  of  significant 
non-zero  singular  values. 

Equations  (2-74)  and  (2-75)  can  be  used  instead  to  obtain  the 
AR  model  matrix  parameters  and  to  determine  the  model  order 
following  the  same  approach  outlined  above.  However,  the  normal 
matrix  in  Equation  (2-74)  is  block  Toeplitz,  and  a  wider  set  of 
alternative  algorithms  can  be  used  to  solve  for  the  matrix 
parameters  (Marple,  1987).  These  two  alternative  systems  of 
equations  (Equations  (2-68)  and  (2-69)  vs.  Equations  (2-74)  and 
(2-75))  can  provide  distinct  results  from  numerical  and  algebraic 
viewpoints.  This  issue  should  be  evaluated  during  Phase  II. 

AR  Model  Identification  Using  the  Raahuveer  Cumulant  Formulation . 

Consider  now  Equation  (2-84),  Raghuveer ' s  linear  system  of 
equations  for  the  AR(L)  model  parameters  {Ak}-  As  in  the  preceding 

case,  Equation  (2-84)  can  be  viewed  as  a  collection  of  J  linear 
systems  of  equations  (corresponding  to  the  J  columns  of  fl) .  In 
this  case,  however,  each  system  consists  of  JLxJL  equations  in  JL 
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unknowns,  and  the  JLxJL  matrix  Cx:L  L  has  full  rank  and  block 
Toeplitz  structure.  The  linear  system  of  Equations  (2-84)  can  be 
solved  using  any  applicable  Toeplitz  equation  solver  algorithm 
(for  example,  Michels  [1991];  Marple  [1987]).  Alternatively, 
Equation  (2-84)  can  be  solved  via  the  pseudoinverse  operator 
implemented  using  the  SVD.  The  SVD  approach  is  preferred  in  cases 
where  model  order  is  unknown  because  the  model  order  determination 
procedure  outlined  for  the  AR  Standard  Cumulant  Formulation  (which 
is  based  on  the  SVD)  can  be  applied  in  this  case  also. 

4 . 4  ARMA  Parameter  Identification 

One  approach  to  multichannel  ARMA  parameter  estimation  is  the 
residual  time  series  method.  In  this  approximate  method  the  AR 
model  order  and  matrix  parameters  are  identified  first,  and  the 
whitening  filter  corresponding  to  the  AR  model  is  applied  to  the 
multichannel  baseband  sequence.  In  theory,  the  whitening  filter 
output  (or  residual)  is  an  MA  sequence.  The  residual  sequence  is 
processed  further  to  estimate  the  MA  model  order  and  matrix 
parameters.  And  the  whitening  filter  corresponding  to  the  MA 
model  is  applied  to  the  residual  output  of  the  first  whitening 
filter.  Also  in  theory,  the  residual  of  the  second  whitening 
filter  is  a  white  sequence.  AR  model  order  and  matrix  parameter 
estimation  is  accomplished  as  described  in  Section  4.3,  except 
that  Equations  (2-89)  and  (2-95)  are  used  instead  of  Equations  (2- 
68)  and  (2-84),  respectively.  MA  model  order  and  matrix  parameter 
estimation  is  accomplished  as  described  in  Section  4.2. 

The  residual  time  series  method  involves  a  simple  concept. 
However,  it  requires  significant  computations  since  two  distinct 
identification  problems  are  solved  in  sequence.  SSC  has  simulated 
the  residual  time  series  method,  and  has  carried  out  various 
simulation-based  analyses.  In  most  cases  this  method  does  not 
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produce  good  results.  Specifically,  SSC  has  observed  that  the  AR 
portion  of  the  procedure  generates  good  matrix  parameter 
estimates,  while  the  MA  portion  generates  poor  estimates  (Section 
5.0).  This  agrees  with  the  observed  performance  of  the  one-slice 
MA  algorithm  applied  to  MA-generated  sequences. 

SSC  has  considered  also  the  hiaher-order  AR  approximation 
method .  which  has  been  used  most  recently  by  Michels  (1991)  for 
complex-valued  vector  sequences.  It  consists  of  fitting  a  higher- 
order  AR  model  to  an  ARMA  sequence.  Preliminary  simulation 
results  suggest  that  the  higher-order  AR  approximation  method 
generates  better  estimates  than  the  residual  time  series  method 
(see  Section  5.1). 

Swami,  Giannakis,  and  Mendel  (1989)  have  proposed  the 
multichannel  M-slice  algorithm  for  ARMA  matrix  parameter 
identification.  This  algorithm  generates  an  estimate  of  the 
matrix  parameters  using  linear  algebra,  without  generating  a 
residual  sequence.  The  algorithm  generates  a  segment  of  the 
impulse  response  sequence  as  an  intermediate  result. 

For  an  ARMA  process  as  defined  in  Equation  (2-25),  the 
imoulse  response  matrix  secruence  {G(n)}  is  defined  as 


(4-20)  y(n)  =  £  G(n-k)u(k)  =  G(k)u(n-k) 

k=-oo  k=-°° 


These  impulse  response  parameters  are  related  to  the  MA  matrix 
parameters  according  to 

k 

(4-21)  BkH  =  X  A" G(k-s)  k  =  1, 2, . . . ,  M 

s=0 
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where  B0  =  lj  and  A0  =  lj  (recall  that  L>M,  and  that  B0  has  full 
rank).  It  follows  that  G(0)  =  lj  . 

The  multichannel  M-slice  algorithm  proceeds  in  two  steps. 
First  the  AR(L)  matrix  parameters,  {Ak},  and  the  first  M+1  impulse 

response  matrices,  {G(n)  I  n  =  0,  1, _ M},  are  estimated  using  a 

generalized  form  of  Equation  (2-89)  .  In  this  first  step  M  1-D 
slices  of  the  output  cumulants  are  used,  thus  giving  rise  to  the 
algorithm's  name.  Next  the  MA  matrix  parameters  are  obtained  using 
Equation  (4-21)  .  This  algorithm  will  be  evaluated  in  Phase  II. 
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5 . 0  SIMULATION-BASED  ANALYSES 


In  Phase  I  SSC  generated  software  routines  to  carry  out 
selected  analyses  and  validate  important  aspects  of  the  HOS-based 
multichannel  detection  methodology  formulated  in  the  program.  The 
software  was  programmed  in  an  Apple  II  Macintosh  processor  using 
FORTRAN  77.  Several  routines,  including  SVD  and  other  matrix 
operations  for  complex-valued  matrices,  were  adapted  from  the 
LINPACK  software  package  (Dongarra  et  al . ,  1979).  The  software 

was  generated  to  carry  out  the  analyses  listed  next. 

(A)  Non-Gauss ian  sequence  generation 

(B)  Identification  algorithm  analyses 

1.  One-slice  MA  identification  algorithm 

2.  Raghuveer  AR  identification  algorithm 

3.  Residual  time  series  ARMA  identification  method 

(C)  Pseudo-innovations  sequence  analyses 

Results  and  conclusions  derived  from  the  analyses  are  summarized 
herein. 

5 . 1  Non-Gaussian  Sequence  Generation 

Two  methods  were  considered  for  generating  colored  (spatial 
and  temporal)  non-Gaussian  vector  sequences.  The  first  method  is 
referred  to  herein  as  the  nonlinear-linear  method.  In  this  method 
a  Gaussian  white  noise  (GWN)  vector  sequence  (zero  spatial  and 
temporal  correlation)  is  generated  and  then  a  zero-memory 
nonlinearity  is  applied  to  each  element  of  the  GWN  vector  sequence 
to  obtain  a  non-Gaussian  white  noise  (NGWN)  vector  sequence.  Then 
the  NGWN  vector  sequence  is  colored,  in  space  and  time,  by 
filtering  with  a  time  series  model.  The  result  is  a  non-Gaussian 
colored  process  (NGCP) .  This  method  produces  a  time  series  with 
known  true  correlation  (spectrum)  and  true  third-order  cumulants, 
but  the  true  PDF  of  the  sequence  is  unknown.  The  linear  filtering 
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process  correlates  the  sequence  and  alters  the  PDF  in  a  manner 
such  that  the  univariate  PDF  of  each  element  of  the  NGCP  is  more 
symmetric,  in  general,  than  the  univariate  PDF  of  the  NGWN . 
Furthermore,  the  univariate  PDF  of  each  element  of  the  NGCP  tends 
to  become  Gaussian  as  the  filter  order  increases.  This  method  is 
useful  for  identification  algorithm  evaluation  when  used  with  low- 
order  shaping  filters.  As  the  filter  order  increases,  longer- 
duration  sequences  are  required  in  order  to  have  reliable 
estimates  of  the  third-order  cumulants  to  use  as  inputs  to  the 
parameter  identification  algorithms.  This  is  due  to  the  fact  that 
third-order  cumulants  of  approximately  symmetric  (or  approximately 
Gaussian)  PDFs  have  small  magnitudes. 

The  second  method  is  referred  to  as  the  1 inear -non 1 inear 
method.  In  this  method  a  GWN  vector  sequence  is  generated  and 
then  colored,  in  space  and  time,  by  filtering  with  a  time  series 
model  to  obtain  a  Gaussian  colored  process  (GCP) .  Then  a  zero- 
memory  nonlinearity  is  applied  to  each  element  of  the  GCP  vector 
to  obtain  a  NGCP.  This  method  produces  a  time  series  with  a  known 
true  univariate  PDF,  but  the  true  correlation  (spectrum)  and  the 
true  third-order  cumulants  of  the  sequence  are  unknown.  The  zero- 
memory  nonlinearity  distorts  the  spatial  and  temporal  statistics 
(of  all  orders)  of  the  sequence.  As  a  result,  the  "true"  linear 
model  that  represents  the  signal  is  unknown.  This  method  is 
useful  for  evaluation  of  cumulant  estimation  techniques.  The 
method  is  useful  also  for  whitening  filter  performance  evaluation, 
specially  for  short-duration  sequences  because  the  sequence  can  be 
generated  to  have  large  third-order  moments  that  are  estimated 
more  accurately.  Also,  some  zero-memory  nonlinear  transformations 
induce  tolerable  distortion  in  the  univariate  second-order 
statistics  (correlation  sequence)  of  the  vector  process.  This  is 
important  for  detection-related  analyses  because  the  detection 
rule  used  in  this  program  is  based  on  second-order  statistics. 
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The  logarithm  is  a  zero-memory  nonlinear  transformation  with 
the  feature  mentioned  above.  Figure  5-1  presents  the  power 
spectrum  of  a  scalar  sequence  generated  by  filtering  a  GWN 
sequence  with  an  AR(2)  .  The  skewness  coefficient  for  the  GCP  at 
the  filter  output  is  0.015,  indicating  a  Gaussian  process.  Figure 
5-2  presents  the  power  spectrum  of  the  sequence  after  application 
of  the  logarithm  operator.  A  comparison  of  the  two  spectra 
indicates  that  the  power  distribution  varies  little  over  the 
frequency  range  (these  figures  represent  single-realization 
spectra);  also,  Figure  5-2  shows  that  the  log  operator  did  not 
introduce  high-frequency  components.  The  skewness  coefficient  of 
the  NGCP  with  spectrum  as  shown  in  Figure  5-2  is  -1.039, 
indicating  its  non-Gaussian  character. 

Other  results  obtained  by  SSC  support  the  utility  of  the  log 
transformation  in  the  context  of  HOS .  In  a  relevant  set  of 
analyses,  a  NGWN  sequence  was  used  to  drive  a  two-channel 
ARMA ( 3 , 2 )  model,  and  GWN  was  added  to  the  ARMA  system  output 
(after  decay  of  initial  transients) .  Raghuveer ' s  AR  algorithm  was 
used  to  identify  the  matrix  parameters  of  AR ( 8 ) ,  AR(10),  .  .  .  , 
and  AR(18)  models.  An  associated  inverse  filter  was  defined  for 
each  of  the  identified  AR  models,  and  the  data  was  processed  with 
each  inverse  filter.  None  of  the  inverse  filters  generated  a 
white  sequence.  This  is  due  to  symmetry  in  the  PDF  of  the 
ARMA (3, 2)  output  induced  by  the  ARMA  transformation  itself.  As  a 
separate  procedure,  the  log  transformation  described  in  Section 
3.3  was  applied  to  the  noisy  ARMA (3, 2)  sequence,  and  Raghuveer ' s 
AR  algorithm  was  applied  to  the  transformed  sequence  to  identify 
the  parameters  of  an  AR(6)  model.  The  inverse  filter  associated 
with  the  AR ( 6 )  model  was  used  to  process  the  transformed  data,  and 
generated  a  white  residual  sequence.  Further  details  on  the 
analyses  and  representative  results  are  presented  next. 
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Figure  5-3  presents  the  procedure  followed  in  the  logarithm 
pre-processing  option  simulation-based  analyses.  The  NGWN  process 
used  as  input  to  the  Data  Generation  step  is  a  zero-mean  vector 
(two-channel)  sequence  of  independent  variables  sampled  from  a 
one-degree-of-freedom  chi-squared  distribution  (the  sample  mean  is 
estimated  and  removed  from  the  chi-squared  sequence) .  This 
sequence  has  a  theoretical  skewness  coefficient  of  21'5  =  2.83. 
The  sample  skewness  coefficient  values  for  the  realization 
selected  for  presentation  herein  are  2.75  and  3.01  for  channels  1 
and  2,  respectively.  A  2,048-point  segment  is  selected  from  the 
steady-state  output  {y(n)}  of  the  ARMA(3,2)  system  (the  true  system 
parameters  are  defined  below) .  Analysis  runs  were  made  also  with 
1,024  data  points,  and  the  results  are  similar  to  those  presented 
herein.  The  sample  skewness  coefficient  values  for  the  {y(n)} 
realization  selected  for  presentation  herein  are  -0.024  and  0.03 
for  channels  1  and  2,  respectively.  These  values  indicate  the 
loss  of  asymmetry  induced  by  the  linear  system.  A  2,048-point 
zero-mean  GWN  vector  sequence  is  added  to  the  steady-state 
ARMA (3,2)  output.  The  noise  variance  in  each  channel  is  selected 
to  result  in  0  dB  SNR. 

Model  Identification  is  carried  out  via  the  higher-order  AR 
approximation  method  using  Raghuveer ' s  AR  algorithm.  Furthermore, 
this  step  is  carried  out  with  and  without  the  logarithm  pre¬ 
processing  option,  as  indicated  in  Figure  5-3.  The  results 
presented  herein  correspond  to  the  case  where  a  two-channel  AR(18) 
model  is  selected  for  the  case  where  pre-processing  is  not  used. 
In  general,  an  AR(18)  model  is  more  than  sufficient  to  represent 
adequately  and  ARMA ( 3 , 2 )  system.  The  logarithm  function  modifies 
the  distribution  that  describes  the  transformed  data.  The  sample 
skewness  coefficient  values  for  the  transformed  data  segment  for 
the  realization  selected  for  presentation  herein  are  -2.05  and 
-1.86  for  channels  1  and  2,  respectively.  These  values  are  close 
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to  the  values  that  describe  the  input  data.  A  two-channel  AR ( 6 ) 
model  is  selected  for  the  case  with  logarithm  pre-processing. 

Residual  Generation  is  accomplished  in  both  cases  using  the 
corresponding  inverse  filter.  Both  inverse  filters  are  two- 
channel  MA  systems  (see  Section  3.5.2) . 

The  true  ARMA(3,2)  parameters  are  (with  A0  =  I2  and  B0  =  I2 )  : 
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This  system  has  six  poles  (number  of  poles  =  number  of  channels  x 
AR  order)  and  six  finite  multivariable  system  zeros  (Rosenbrock, 
1970).  The  system  poles  are  located  at:  {0.289, -0.436, -0.552  ±  j  0.778, 
0.218  ±  j  0.549},  and  the  system  zeros  are  located  at:  {0.7,  0.7,  0.7,  0.7, 
0.0,  0.0}.  This  system  is  minimum  phase  since  all  the  poles  and 
zeros  are  inside  the  unit  circle.  Thus,  the  true  inverse  system 
(whitening  filter)  is  stable  also. 

Residual  auto-correlation  sequences  for  a  single  realization 
of  the  channel  1  output  for  each  of  the  two  cases  are  presented  in 
Figures  5-4  and  5-5.  Channel  2  residual  auto-correlations  are 
similar.  In  each  figure  the  complete  4,096-point  auto-correlation 
sequence  is  shown.  The  central  segment  of  the  auto-correlation  is 
shown  also  in  order  to  highlight  the  correlation  detail.  Notice 
that  the  residual  is  colored  for  the  case  without  pre-processing 
(Figure  5-4),  and  white  for  the  case  with  logarithm  pre-processing 
(Figure  5-5)  .  These  figures  are  representative  of  all  the 
realizations  generated  in  the  analysis. 
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Figure  5-3 .  Set-up  for  logarithmic  pre-processing  option 

simulation-based  analyses. 
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amplitude 


Auto-correlation  of  channel  1  residual  of  ARMA(3,2)  process 
filtered  with  whitening  filter  for  AR(1 8)  model  (SNR  =  0  dB) 


(a)  Complete  auto-correlation  sequence  for  Channel  1  residual. 


Auto-correlation  of  channel  1  residual  of  ARMA(3,2)  process 
filtered  with  whitening  filter  for  AR(1 8)  model  (SNR  =  0  dB) 


(b)  Central  segment  of  auto-correlation  sequence  for  channel  1 

residual . 

Figure  5-4.  Auto-correlation  of  residual  for  AR(18)  model 
(without  pre-processing)  of  ARMA(3,2)  plus  GWN  data. 


amplitude 


Auto-correlation  of  channel  1  residual  of  log  of  ARMA(3,2)  process 
filtered  with  whitening  filter  for  AR(6)  model  (SNR  =  0  dB) 


(a)  Complete  auto-correlation  sequence  for  Channel  1  residual. 


Auto-correlation  of  channel  1  residual  of  log  of  ARMA(3,2)  process 
filtered  with  whitening  filter  for  AR(6)  model  (SNR  =  0  dB) 


lag  index,  m 

(b)  Central  segment  of  auto-correlation  sequence  for  channel  1 

residual . 

Figure  5-5.  Auto-correlation  of  residual  for  AR(6)  model  (with 
log  pre-processing)  of  ARMA(3,2)  plus  GWN  data. 
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5 . 2  Identification  Algorithms  Analyses 


Software-based  analyses  were  carried  out  to  validate 
analytical  derivations  and  to  determine  algorithm  performance. 
The  results  presented  herein  are  representative  of  all  the  results 
obtained  in  Phase  I.  The  additive  noise,  {w(n)},  is  zero  for 
Examples  5-1  through  5-5. 

5.2.1  One-Slice  MA  Identification  Algorithm 

In  most  of  the  test  cases  considered  the  MA  algorithm 
performed  poorly  in  an  absolute  sense,  as  well  as  in  comparison 
with  the  AR  algorithm.  The  MA  algorithm  generates  non-minimum- 
phase  zeros  more  often  than  the  AR  algorithm,  even  when  the  true 
zeros  are  minimum-phase.  This  performance  is  evidenced  in  the 
examples  presented  next  for  a  MA(2)  model  (the  results  can  be 
compared  directly  with  the  results  from  the  AR  examples  presented 
in  Section  5.2.2)  . 

Two  different  methods  are  used  to  obtain  the  matrix  parameter 
estimates  in  the  examples.  In  the  first  method  (Method  I),  the  MA 
parameter  estimates  are  obtained  using  a  single  realization  of  104 
time  samples.  In  the  second  method  (Method  II),  ten  realizations 
of  103  time  samples  each  are  processed  to  obtain  ten  sets  of 
estimated  MA  parameters;  the  average  of  the  ten  sets  of  parameters 
is  calculated  and  presented  herein.  Notice  that  the  total  number 
of  samples  used  is  the  same  in  both  methods . 

With  respect  to  whitening  filter  performance,  the  locations 
of  the  zeros  of  the  transfer  function  of  the  identified  MA  model 
determine  the  system  dynamics.  In  these  examples,  the  values 
listed  as  "zeros"  are  the  roots  of  the  polynomial  (3(z),  as  defined 
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in  Equation  (2-41) .  Also,  the  normalized  error  in  the  estimate  of 
the  ith  zero,  denoted  as  §Zj,  is  defined  as 


:5-i: 


8Z;  =  |Zi'Z|1 


Z; 


where  Zj  denotes  the  true  ith  zero  and  Zj  denotes  its  estimate.  This 
error  is  converted  to  a  percent  error  by  multiplication  by  100%. 


Examole  5- 

-_1.  The  true  MAC 

b"  = 

[-1.0  0 

0 

1 

b 

B2  = 

0.25  0 

0  0.25. 

Method  I  Estimates 

(single 

Bi  = 

-0.9991 

0.0078 

0.0073 

-1.0383  - 

= 

0.2425 

-0.0071  1 

-0.0080 

0.2561  J 

zeros:  {0.5,  0.5,  0.5,  0.5} 


zeros:  {0.388,  0.435,  0.577,  0.637} 


The  errors  in  the  location  of  the  four  real-valued  zeros  in  the 
order  listed  above  are:  22.4%,  13.0%,  15.4%,  and  27.4%. 


Method  II 


Estimates 


(average  of  ten  realizations) : 


bH  _  [  -0.9968  0.0111 

1  L  -0.0007  -0.8993 


zeros:  {0.445  ±  j  0.1 54,  0.503  ±  j  0.232} 


0.3036  0.0094 

-0.0062  0.2239 
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The  error  in  the  location  of  the  zeros  is  32.7%  for  the  first 
complex  conjugate  pair,  and  46.4%  for  the  second  complex  conjugate 
pair.  Some  of  the  zeros  of  the  matrix  parameters  estimated  for 
the  individual  realizations  have  a  larger  error;  the  maximum  error 
is  91.8%  for  the  complex  conjugate  pair  estimated  at  0.48  ±  j 
0.46.  These  large  errors  are  evident  in  Figure  5-6,  a  scatter 
plot  of  the  forty  zeros  estimated  in  the  ten  realizations .  All 
estimated  sets  of  zeros  are  plotted  in  a  single  set  of  axes  since 
all  the  true  zeros  are  at  the  same  location,  Z  =  0.5.  Error 
statistics  for  the  estimated  zeros  are  presented  in  Tables  5-1  and 
5-2.  All  estimated  zeros  are  inside  the  unit  circle,  as  evidenced 
in  Figure  5-6. 


ZERO  LOCATIONS  (MA(2)  SYSTEM) 


real 


Legend: 

True  zeros:  * 
Estimated  zeros:  o 


Figure  5-6.  Scatter  plot  for  the  zeros  of  all  ten  realizations  in 

Example  5-1  (Method  II) . 
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Example  5-2.  The  true  MA(2)  parameters  are  (with  B0  =  12)  : 

gH  _  -1.8  0 

1  “l  0  -1.8  J  zeros:  {0.9,  0.9,  0.9,  0.9} 

gH  _  0.81  0 

2  “  0  0.81 

Method  I  Estimates  (single  realization) : 

bh  _  r  -1.8100  0.0222 

1  1-0.0035  -1.8237  J  zeros:  {0.856,  0.909  ±  j  0.145,  0.969} 

bh  _  [  0.8313  -0.0052 

2  -0.0108  0.8457  . 

The  error  in  the  locations  of  the  zeros  is  4.9%  for  the  first 
real-valued  zero,  16.1%  for  the  complex  conjugate  pair,  and  7.7% 
for  the  second  real-valued  zero. 

Method  II  Estimates  (average  of  ten  realizations) : 

gH  =  [  -1.7070  -0.0047' 

1  i  0.0243  -1.6572  J  zeros:  {0.826  ±  j  0.260,  0.856  ±  j  0.304} 

gH  =  [  0.8209  -0.0003 

2  -0.0148  0.7542  - 

The  error  in  the  locations  of  the  zeros  is  30.0%  for  the  first 
complex  conjugate  pair,  and  34.1%  for  the  second  complex  conjugate 
pair.  The  zeros  of  several  of  the  matrix  parameters  estimated  for 
the  individual  realizations  have  a  larger  error,  and  eleven  (out 
of  forty)  zeros  are  non-minimum-phase.  This  is  evidenced  in 
Figure  5-7,  which  is  a  scatter  plot  of  the  forty  zeros  estimated 
in  the  ten  realizations.  The  largest  error  is  56.7%,  and 
corresponds  to  a  non-minimum-phase  zero  at  1.41.  Figure  5-7  is  in 
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the  same  scale  as  Figure  5-6.  Notice  that  the  dispersion  of  zeros 
is  about  the  same  in  both  figures . 


ZERO  LOCATIONS  (MA(2)  SYSTEM) 


Figure  5-7.  Scatter  plot  for  the  zeros  of  all  ten  realizations  in 

Example  5-2  (Method  II) . 


Example  5-3 .  The  true  MA(2)  parameters  are  (with  B0  =  I2 )  : 


B 


H 

1 


-2.2  0 
0  -2.2 . 


1.21  0 
0  1.21. 


zeros: {1.1,  1.1,  1.1,  1.1} 


Notice  that  this  MA  system  is  non-minimum  phase  (HOS-based  methods 
can  identify  non-minimum  phase  systems,  whereas  covariance-based 
methods  cannot).  As  mentioned  earlier,  non-minimum  phase  systems 


result  in  unstable  innovations  filters.  This  example  is  presented 
to  demonstrate  the  identification  capabilities  of  HOS-based 
methods  in  a  more  general  context. 


Method  I  Estimates  (single  realization) : 


B 


H 

1 


-2.2411  0.0297 

-0.0119  -2.2144  . 


zeros:  {1.027,  1.110  ±j  0.173,  1.21} 


B 


H 

2 


1.2535  -0.0089 

0.0215  1.2493  . 


The  error  in  the  locations  of  the  zeros  is  6.6%  for  the  first 
real-valued  zero,  15.8%  for  the  complex  conjugate  pair,  and  10.0% 
for  the  second  real-valued  zero.  All  estimated  zeros  are  non¬ 
minimum  phase. 


Method  II  Estimates  (average  of  ten  realizations) : 


B 


H 

1 


-2.0595  -0.0213 

0.0409  -2.0469  . 


zeros:  {1 .006  ±  j  0.31 6,  1 .047  ±  j  0.356} 


B 


H 

2 


1.1852  0.0079 

-0.0244  1.1480  . 


The  error  in  the  location  of  the  zeros  is  30.0%  for  the  first 
complex  conjugate  pair,  and  32.7%  for  the  second  complex  conjugate 
pair.  The  zeros  of  various  parameter  estimates  for  the  individual 
realizations  have  a  larger  error;  the  largest  error  is  67.1%,  and 
corresponds  to  a  real-valued  zero  at  1.84.  This  is  evidenced  in 
Figure  5-8,  which  is  a  scatter  plot  of  the  forty  zeros  estimated 
in  the  ten  realizations.  Note  that  out  of  the  forty  estimated 
zeros,  eleven  are  estimated  incorrectly  to  be  inside  the  unit 
circle . 
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real 

Figure  5-8.  Scatter  plot  for  the  zeros  of  all  ten  realizations  in 

Example  5-3  (Method  II) . 

Comments .  The  statistics  for  the  zeros  of  the  parameter  estimates 
of  the  ten  realizations  in  Method  II  for  the  three  examples  are 
listed  in  Table  5-1.  These  results  indicate  that  "slow"  zeros 
(close  to  the  unit  circle)  are  estimated  more  accurately  than 
"fast"  zeros.  This  feature  is  common  to  all  identification 
algorithms  (Roman  and  Davis,  1993) .  Figures  5-6  through  5-8 
suggest  that  the  magnitude  of  the  estimated  zeros  is  biased 
towards  a  magnification  error. 

Table  5-2  lists  the  errors  in  the  zeros  of  the  parameters 
estimated  via  Methods  I  and  II  for  the  three  MA(2)  examples.  In 
all  three  examples  the  single  realization  method  generated  better 
zero  estimates  than  the  averaging  method.  This  observation  admits 
the  following  explanation.  In  the  single  realization  method  all 
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the  averaging  takes  place  in  the  estimation  of  the  cumulants, 
prior  to  processing  with  the  parameter  identification  algorithm. 
In  the  averaging  method  part  of  the  averaging  takes  place  after 
processing  with  the  parameter  identification  algorithm.  The 
identification  algorithm  propagates  the  input  errors,  and  the 
post-processing  averaging  is  less  effective.  Thus,  for  the  MA 
algorithm,  it  is  preferable  to  use  more  accurate  cumulant 
estimates  than  to  use  post-algorithm  averaging. 


TRUE 

ZERO 

(REPEATED) 

mm 

mm 

ERROR 

MEAN 

(%) 

ERROR 
STD.  DEV. 

<%) 

ERROR 

MEDIAN 

.(%) 

EX. 

5-1 

in 

o 

4.7 

91.8 

46.9 

21.9 

48.3 

EX. 

5-2 

0.9 

mm 

56.7 

24.0 

1  1.9 

23.5 

EX. 

5-3 

1.1 

9.3 

67.1 

23.5 

1  1.6 

19.7 

Table  5-1.  Error  statistics  for  the  zeros  of  all  ten  realizations 
in  the  three  MA  examples  (Method  II) . 


ERRORS  USING  METHOD  1 

(%) 

ERRORS  USING  METHOD  II 

(%) 

EXAMPLE  5-1 

22.4,  13.0,  15.4,  27.4 

32.7,  32.7,  46.4,  46.4 

EXAMPLE  5-2 

4.9,  16.1,  16.1,  7.7 

30.0,  30.0,  34.1,  34.1 

EXAMPLE  5-3 

6.6,  15.8,  15.8,  10.0 

30.0,  30.0,  32.7,  32.7 

Table  5-2.  Zero  estimate  errors  for  the  three  MA(2)  examples 

using  both  methods. 
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5.2.2  Raghuveer  AR  Identification  Algorithm 

The  AR  algorithm  generates  non-minimum-phase  poles  less  often 
than  the  MA  algorithm  (the  values  referred  to  as  "poles"  are  the 
roots  of  the  polynomial  a(z),  as  defined  in  Equation  (2-40))  .  This 

is  evidenced  in  Examples  5-4  and  5-5  for  an  AR(2)  model  with  true 
poles  located  at  the  same  values  as  the  zeros  in  Examples  5-1  and 
5-2,  respectively.  In  order  to  allow  direct  comparisons,  the  true 
AR ( 2 )  parameters  are  selected  equal  to  the  true  MA(2)  parameters 
in  the  respective  examples.  Also,  the  normalized  error  in  the 
estimate  of  the  ith  pole,  denoted  as  5p,,  is  defined  analogous  to 
5z,  in  Equation  (5-1). 

Example  5-4.  The  true  AR(2)  parameters  are  the  same  as  the  true 
MA(2)  parameters  of  Example  5-1;  that  is,  A2  =  B2,  A1  =  B1 ,  and  A0 
=  l2.  Also,  the  true  poles  are  located  at  {0.5,  0.5,  0.5,  0.5} . 

Method  I  Estimates  (single  realization) : 

aH  _  [  -1.0065  0.0060 

1  L  0.0026  -1.0029  J  poles:  {0.501  ±j  0.048,  0.504  ±j  0.066} 

ah  _  [  0.2586  -0.0001 

2  i  -0.0023  0.2528  . 

The  error  in  the  first  complex  conjugate  pole  pair  is  9.6%,  and 
the  error  in  the  second  complex  conjugate  pole  pair  is  13.1%. 

Method  II  Estimates  (average  of  ten  realizations): 

ah  _  [  -0.9108  -0.0048  ' 

1  L  -0.0051  -0.9147  J  poles:  {0.454  ±j  0.164,  0.459  ±j  0.145} 

ah  =  [  0.2313  -0.0023  ' 

2  l  0.0005  0.2332  . 


109 


The  error  in  the  first  complex  conjugate  pole  pair  is  34.1%,  and 
the  error  in  the  second  complex  conjugate  pole  pair  is  30.1%. 
These  errors  are  slightly  larger  than  the  mean  error  for  all 
forty  matrix  parameters  estimated  for  the  individual 
realizations,  which  is  21.7%  for  the  complex  conjugate  pole  pair 
estimated  at  0.49  ±  j  0.53  (see  Table  5-3) .  From  Table  5-3,  the 
maximum  error  is  44.8%,  and  corresponds  to  the  pole  estimated  at 
0.28.  A  scatter  plot  of  the  forty  poles  estimated  for  the  ten 
realizations  is  presented  in  Figure  5-9.  Notice  that  the 
dispersion  in  Figure  5-9  is  significantly  less  than  the 
dispersion  in  Figure  5-6.  Notice  also  that  the  complex-valued 
pole  estimates  appear  along  a  circle  of  magnitude  0.5. 


POLE  LOCATIONS  (AR(2)  SYSTEM) 


Legend: 

True  poles:  * 
Estimated  poles:  x 


Figure  5-9.  Scatter  plot  for  the  poles  of  all  ten  realizations 

in  Example  5-4  (Method  II) . 
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Example  5-5 .  The  true  AR(2)  parameters  are  the  same  as  the  true 
MA  ( 2 )  parameters  of  Example  5-2;  that  is,  A2  =  B2,  A-|  =  B-| ,  and  A0 
=  l2.  Also,  the  true  poles  are  located  at  {0.9,  0.9,  0.9,  0.9}. 


Method  I  Estimates  (single  realization) : 


-1.7975  0.0070 

0.0121  -1.8044 


poles:  {0.870,  0.921 ,  0.906  ±  j  0.038} 


0.8075  -0.0071 

-0.0119  0.8151 


The  two  real-valued  poles  have  an  error  of  3.3%  and  2.3%, 
respectively,  and  the  error  in  the  complex  conjugate  pole  pair 
is  4.3%. 


Method  II  Estimates  (average  of  ten  realizations) : 


A 


H 

1 


-1.6350  -0.0248 

0.0099  -1.6432  . 


poles:  {0.813  ±  j  0.269,  0.826  ±  j  0.253} 


0.7355  0.0242 

-0.0108  0.7439 


The  error  in  the  first  complex  conjugate  pole  pair  is  31.4%,  and 
the  error  in  the  second  complex  conjugate  pole  pair  is  29.3%. 
Figure  5-10  is  a  scatter  plot  of  the  forty  poles  estimated  for  the 
ten  realizations.  The  poles  of  the  estimated  individual  matrix 
parameters  have  smaller  error  than  the  poles  of  the  averaged 
matrix  parameters  given  above.  In  fact,  the  maximum  error  is 
13.5%  for  the  pole  estimated  at  0.78  (see  Table  5-3)  .  This  is 
possibly  due  to  the  fact  that  the  poles  are  a  nonlinear  function 
of  the  averaged  matrix  parameters.  From  Figure  5-10,  the  poles  of 
all  ten  sets  of  estimated  matrix  parameters  are  minimum-phase. 


Ill 


And  the  dispersion  in  Figure  5-10  is  much  less  than  in  Figures  5-7 
(MA  case)  and  5-9.  Additionally,  in  this  example  the  dispersion 
does  not  exhibit  a  pattern. 


POLE  LOCATIONS  (AR(2)  SYSTEM) 


real 

Figure  5-10.  Scatter  plot  for  the  poles  of  all  ten  realizations 

in  Example  5-5  (Method  II) . 

Comments .  The  statistics  for  the  poles  of  the  parameter  estimates 
of  the  ten  realizations  in  Method  II  for  Examples  5-4  and  5-5  are 
listed  in  Table  5-3.  These  results  also  indicate  that  "slow" 
poles  (close  to  the  unit  circle)  are  estimated  more  accurately 
than  "fast"  poles,  as  is  the  case  for  the  zeros.  The  magnitude 
error  bias  observed  in  the  MA  cases  is  present  also  in  these  AR 
cases  using  Raghuveer '  s  algorithm,  but  to  a  significantly  lesser 
degree . 
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Table  5-4  lists  the  errors  in  the  poles  of  the  parameters 
estimated  via  Methods  I  and  II  for  the  two  AR(2)  examples.  In 
both  examples  the  single  realization  method  generated  better  pole 
estimates  than  the  averaging  method,  which  admits  the  same 
explanation  as  in  the  MA  case.  Thus,  for  the  AR  algorithm  it  is 
preferable  also  to  use  more  accurate  cumulant  estimates  than  to 
use  post-algorithm  averaging. 


BM 

MINIMUM 

ERROR 

(%> 

MAXIMUM 

ERROR 

(%) 

ERROR 

MEAN 

(%) 

ERROR 
STD.  DEV. 

(%) 

ERROR 

MEDIAN 

(%) 

EX.  5-4 

0.5 

10.8 

44.8 

21.7 

7.7 

22.5 

EX.  5-5 

0.9 

2.2 

13.5 

7.1 

2.6 

6.7 

Table  5-3.  Error  statistics  for  the  poles  of  all  ten  realizations 

in  the  two  AR  examples  (Method  II) . 


ERRORS  USING  METHOD  1 

(%) 

ERRORS  USING  METHOD  II 

(%) 

EXAMPLE  5-4 

9.6,  9.6,  13.1,  13.1 

34.1,  34.1,  30.1,  30.1 

EXAMPLE  5-5 

3.3,  2.3,  4.3,  4.3 

31.4,  31.4,  29.3,  29.3 

Table  5-4.  Pole  estimate  errors  for  the  two  AR(2)  examples  using 

both  methods . 


Direct  comparison  of  Figures  5-6  and  5-7  with  Figures  5-9  and 
5-10,  and  of  Tables  5-1  and  5-2  with  Tables  5-3  and  5-4  indicates 
that  the  AR  algorithm  generates  much  better  estimates  than  the  MA 
algorithm  considering  all  relevant  performance  measures. 
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Effect  of  Additive  Gaussian  White  Noise  on  AR  Estimates .  The 


effect  of  additive  Gaussian  noise  is  considered  in  Example  5-6. 
The  example  is  for  a  specific  set  of  conditions  using  Raghuveer 1 s 
AR  algorithm,  but  the  results  are  representative  of  parameter 
identification  performance  in  noise  for  AR  and  MA  algorithms. 

The  normalized  error  defined  via  Equation  (5-1)  is  useful  in 
the  determination  of  dynamic  system  response  effects  .  An 
alternative  error  definition  is  required  to  assess  parameter 
estimation  performance.  A  useful  error  measure  is  defined  as 
follows.  Let  Ci;k  denote  the  ith  column  of  Ak,  and  let  C:.k  denote 

its  estimate  obtained  using  Raghuveer ' s  AR  algorithm.  The 
normalized  average  column  error,  denoted  as  5c,  is  defined  as 

<5-2.  &  = 

JL  i=i  k=1  I  — i;k 

This  error  is  converted  to  a  percent  error  by  multiplication  by 
100%.  In  the  example,  be  is  computed  for  each  realization  as  a 

function  of  noise  level. 

Example  5 - 6  .  The  signal  is  a  two-channel  AR(2)  process  with 
(four)  repeated  poles  at  0.5  and  zero  coupling  between  channels, 
as  in  Example  5-4.  Gaussian  white  noise  is  added  to  the  AR(2) 
process.  Sequence  duration  is  5xl04  time  samples  for  each 
realization.  Raghuveer 1 s  AR  parameter  estimation  algorithm  is 
used  to  identify  the  AR  matrix  parameters.  The  normalized  average 
column  percent  error  is  plotted  in  Figure  5-11.  Notice  that  for 
SNR  >  7  dB  the  identification  error  is  close  to  its  noise-free 
value  (the  error  is  less  than  1%  off  the  noise-free  asymptote) . 
This  is  a  consequence  of  the  insensitivity  of  third-order 
cumulants  to  additive  Gaussian  noise.  The  error  increases 
dramatically  for  SNR  <  7  dB .  The  SNR  value  at  which  the  plot 
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starts  to  deviate  from  the  noise-free  asymptote  is  a  function  of 
sequence  duration  also. 

Normalized  rms  error  for  AR  coefficient  estimates 


SNR  (dB) 


Figure  5-11.  Normalized  average  column  error  for  AR  matrix 
parameters  as  a  function  of  SNR  (additive  Gaussian  noise) . 

5.2.3  Residual  Time  Series  ARMA  Identification  Method 

The  residual  time  series  ARMA  identification  method  did  not 
produce  good  results  in  most  of  the  cases  considered  in  Phase  I. 
Simulation-based  analyses  indicate  that  the  AR  algorithm  generates 
good  matrix  parameter  estimates,  while  the  MA  algorithm  generates 
poor  estimates.  As  evidenced  in  Section  5.2.1,  the  MA  algorithm 
generates  non-minimum-phase  estimates  in  many  cases,  even  in  the 
identification  of  the  parameters  of  a  true  minimum-phase  MA 
process.  In  the  residual  method  the  situation  is  complicated 
further  by  the  fact  that  the  MA  algorithm  has  to  deal  with 
limited-duration  residuals  (from  the  AR  inverse  filter)  that  are 
not  a  true  MA  process  (unless  the  identified  AR  parameters  are  the 
true  parameters) . 
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An  alternative  approach  is  to  fit  a  higher-order  AR  model  to 
the  ARMA  sequence.  This  approach  has  been  considered  by  Michels 
(1991)  for  Case  (2-20d),  the  multichannel  Gaussian  case. 
Preliminary  results  obtained  by  SSC  suggest  that  such  an  approach 
is  preferable  over  the  residual  method  in  the  non-Gaussian  case 
(see,  for  example,  the  results  in  Section  5.1) . 

5 . 3  Pseudo-Innovations  Sequence  Analyses 

The  capability  of  the  HOS-based  detection  methodology  to 
discriminate  between  the  two  hypotheses  was  established  in  Phase 
I.  Simulation  analyses  were  carried  out  for  non-Gaussian  colored 
sequences  in  additive  white  noise  generated  using  the  linear- 
nonlinear  method,  with  the  log  function  as  the  zero-memory 
nonlinearity.  GWN  was  added  to  the  GCP  prior  to  applying  the  log 
to  the  sequence.  The  results  show  that  the  output  of  the  H0  filter 
is  white  when  the  null  hypothesis  is  true,  {x(n)  =  c(n)  +  w(n)} ,  and  the 
output  of  the  H0  filter  is  colored  when  the  alternative  hypothesis 
is  true,  {x(n)  =  s(n)  +  c(n)  +  w(n)} .  Similarly,  the  output  of  the  H1 
filter  is  white  when  the  alternative  hypothesis  is  true,  {x(n)  =  s(n) + 
c(n)+w(n)},  and  the  output  of  the  H-j  filter  is  colored  when  the  null 
hypothesis  is  true,  {x(n)  =  c(n)  +  w(n)} .  Given  these  observations,  an 
appropriate  detection  rule  (such  as  the  approximate  rule 
considered  in  Section  3.4)  can  discriminate  between  the  two 
hypotheses  based  on  the  differences  in  the  auto-correlation  of  the 
pseudo- innovations .  A  representative  sample  of  this  type  of 
result  is  provided  in  Example  5-7.  Detailed  analyses  will  be 
carried  out  in  Phase  II  to  quantify  detection  performance. 

Example  5-7  .  Both  clutter  and  signal  were  generated  using  the 
linear-nonlinear  method.  Specifically,  the  clutter  is  a  two- 
channel  AR(2)  process  driven  by  GWN,  and  the  signal  is  also  a  two- 
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channel  AR ( 2 )  process  driven  by  GWN.  The  true  AR(2)  matrix 
parameters  for  clutter  and  signal  are  (these  are  the  real  part  of 
the  complex-valued  model  parameters  in  [Roman  and  Davis,  1993]): 


CLUTTER 

SIGNAL 

-1.0430  0.0 

ah  - 

Mi  “ 

'1.6290 

0.0 

0.0  -1.0430  . 

.  0.0 

1.6290. 

0.4900  0.0 

Ap  = 

0.8099 

0.0 

.  0.0  0.4900. 

d 

0.0 

0.8099- 

GWN  was  added  to  the  AR  sequences  to  generate  null  and  alternative 
hypotheses  data,  as  appropriate,  at  an  SNR  level  of  0  dB .  The 
addition  of  GWN  converts  the  sequence  into  an  ARMA  sequence.  A 
sequence  duration  of  1,024  time  samples  (after  initial  transients 
decay)  was  selected.  The  Gaussian-distributed  data  (including  the 
additive  noise)  was  transformed  using  the  logarithm  transformation 
as  described  in  Section  3.3.  This  completes  the  data  generation 
step  according  to  the  linear-nonlinear  method  (Section  5.1) .  An 
alternative  procedure  is  to  add  the  GWN  after  the  log  is  applied 
to  the  AR  data.  This  should  result  in  better  perf ormance . 

Raghuveer ' s  AR  identification  algorithm  was  used  to  generate 
AR  matrix  parameter  estimates  for  a  sixth-order  AR  model  for  the 
null  hypothesis  case  (clutter  +  noise) .  Figure  5-12  presents  the 
auto-correlation  function  of  the  pseudo-innovations  e(nlH0),  the 
residual  at  the  output  of  the  sixth-order  Hg  filter  driven  with  H0 
data.  The  residual  is  uncorrelated,  as  expected.  Figure  5-13 
presents  the  auto-correlation  function  of  the  residual  at  the 
output  of  the  sixth-order  H0  filter  driven  with  H-|  data.  As 

expected,  the  residual  in  Figure  5-13  is  correlated  over  several 
lags.  This  difference  is  the  mechanism  by  which  discrimination 
between  the  two  hypotheses  takes  place. 
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amplitude 


Auto-correlation  of  channel  1  residual  for  case  of  null  hypothesis  data 

filtered  with  null  hypothesis  filter  (6th-order,  SNR=0dB) 
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Figure  5-12.  Central  segment  of  pseudo -innovations  auto¬ 
correlation  (H0  data  filtered  with  H0  filter)  for  Example  5-7 
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Auto-correlation  of  channel  1  residual  for  case  of  alternative  hypothesis  data 
filtered  with  null  hypothesis  filter  (6th-order,  SNR=0dB) 
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Figure  5-13  .  Central  segment  of  pseudo- innovations  auto¬ 
correlation  (H-|  data  filtered  with  H0  filter)  for  Example  5-7 


6 . 0  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  program  the  feasibility  of  using  higher-order 
cumulants  for  radar  target  detection  was  established  in  a  general 
context,  as  well  as  in  the  specific  context  of  surveillance  radar 
applications.  With  respect  to  surveillance  radar,  SSC  identified 
several  operational  scenarios  where  non-Gaussian  targets  are 
present  in  a  non-Gaussian  clutter  environment.  Additionally,  at 
least  two  different  candidate  problems  of  interest  to  RL  were 
identified:  radar  arrays  and  pre-detection  radar  fusion.  Emphasis 
was  placed  in  surveillance  radar  arrays,  which  include  the 
space/time  processing  application.  In  this  context,  a  new 
methodology  for  model-based  multichannel  detection  was  formulated 
using  third-order  cumulants  to  identify  the  model  parameters  for 
models  in  the  time  series  class.  This  methodology  builds  on  the 
work  of  Metford  (1984)  and  of  Michels  (1991)  by  the  significant 
inclusion  of  higher-order  statistics. 

The  SSC  methodology  can  be  applied  for  single-cell  detection, 
which  corresponds  to  the  cases  wherein  the  temporal  statistics  of 
the  radar  return  are  non-Gaussian.  The  methodology  can  be  applied 
also  for  global  detection,  wherein  the  spatial  statistics  of  the 
radar  return  are  non-Gaussian.  In  these  applications  the  use  of 
third-order  cumulants  to  identify  model  parameters  provides 
immunity  to  additive  Gaussian  noise  and  interference  sources 
because  all  higher-order  cumulants  are  zero  for  Gaussian- 
distributed  processes. 

Simulation-based  analyses  were  carried  out  to  validate  key 
aspects  of  the  methodology,  and  to  establish  the  feasibility  of  an 
HOS-based  approach  to  model-based  multichannel  detection.  Results 
obtained  to  date  suggest  the  capability  to  discriminate  between 
target -present  and  target-absent  hypotheses  using  the  suboptimal 
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and  approximate  detection  rule  proposed  by  Metford  (1984)  for  non- 
Gaussian  signals. 

In  Phase  I  several  key  technical  issues  were  addressed,  and 
candidate  resolution  approaches  were  identified.  Those  issues  are 
discussed  in  this  report,  and  are  summarized  in  Table  6-1.  Most 
of  those  issues  will  be  addressed  further  during  Phase  II. 

Several  technical  tasks  have  been  identified  also  for  further 
research  and  development  in  Phase  II.  These  tasks  are  summarized 
next . 

Processor  Requirements  Definition.  Determination  of  the  potential 
of  the  SSC  methodology  for  radar  system  applications  requires  the 
establishment  of  requirements  for  the  radar  problems  of  interest 
at  RL .  This  includes  space/time  processing  in  a  surveillance 
radar  array  system,  and  fusion  of  data  from  multiple  distinct 
radar  systems.  Similarly,  requirements  must  be  identified  for 
non-defense  applications  of  interest. 

Additional  Analyses  and  Detailed  Algorithm  Formulation .  The 
analyses  listed  below  are  required  in  order  to  generate  a  detailed 
algorithm  definition  and  to  assess  the  SSC  methodology  in  the 
context  of  the  applications  requirements.  These  analyses  support 
the  technical  issues  listed  in  Table  6-1. 

•  Alternative  matrix  parameter  identification  algorithms  (such 
as  the  multichannel  M-slice  algorithm  for  ARMA  parameter 
estimation)  will  be  analyzed  to  select  the  one  which  best 
satisfies  the  trade-offs  associated  with  performance  and 
other  relevant  criteria. 
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ISSUE 


COMMENTS 


CANDIDATE  RESOLUTION 
APPROACHES 


Possible 

symmetry  of 

marginal  PDF  of 

the  quadrature 

components. 

Not  an  issue  in  the  cases  (if  any) 

where  the  marginal  PDF  of  the 

quadrature  components  is 

asymmetric.  (The  marginal  PDF  is 

unknown  for  most  non-Gaussian 

quadrature  components.) 

A.  Apply  zero-memory  nonlinear 

transformation  to  the  quadrature 

components. 

B.  Use  formulation  based  on  fourth-order 

cumulants. 

Non-minimum 

•  Major  issue  only  for  MA  models. 

A.  Average  estimates  over  several 

phase  model 

•  Impact  on  ARMA  models  needs  to 

resolution  cells. 

parameter 

be  assessed  further. 

B.  Use  AR  (or  possibly  ARMA)  models  only. 

estimates. 

C.  Use  non-causal  formulation. 

Complex-valued 

Restrictions  are  enhanced  for  non- 

A.  Use  PDF  representation  based  on 

probabilistic 

representation  is 

restrictive. 

Gaussian  sequences. 

concatenated  real-valued  vectors. 

Detection  rule  for 

Sub-optimal  approximations  may 

A.  Use  Metford's  approximate  detection 

non-Gaussian 

provide  satisfactory  performance. 

rule  for  non-Gaussian  sequences. 

sequences. 

B.  Use  non-parametric  detection  rule. 

•  Generation  of  temporal 

A.  Apply  zero-memory  nonlinear 

multichannel  processes  with 

transformation  to  the  quadrature 

Non-Gaussian 

specified  multivariate  PDF  is  an 

components. 

multichannel  data 

unsolved  problem. 

B.  Use  SIRP-based  data  generation  for 

generation  for 

•  Linear  filtering  of  a  non-Gaussian 

global  domain  detection  problems. 

performance 

process  induces  Gaussianity  or 

C.  Investigate  alternative  techniques  such 

evaluation. 

symmetry  in  the  PDF.  This 

as  extension  of  the  Liu  and  Munson 

increases  the  size  of  the  data  base 

required  to  obtain  satisfactory 

cumulant  estimates. 

(1982)  approach  to  vector  processes. 

Table  6-1.  List  of  key  technical  issues  and  candidate  resolution 

approaches . 
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•  Giannakis  and  Delopoulos  (1990)  have  demonstrated  that  the 
covariance  sequence  can  be  estimated  from  the  third-order 
cumulants .  Thus,  it  is  possible  to  define  an  identification 
approach  which  combines  the  insensitivity  to  additive 
Gaussian  noise  offered  by  third-order  cumulants  with  the 
robustness  of  covariance-based  model  identification  methods 
for  parametric  models  (Michels,  1991;  Roman  and  Davis,  1993)  . 
This  combined  approach  will  be  investigated,  with  emphasis  on 
the  accuracy  of  the  HOS-based  covariance  sequence  estimate. 

•  Model  order  selection  criteria  for  on-line  and  off-line 
decisions  will  be  evaluated  and  the  preferred  criterion  will 
be  adopted.  Criteria  to  be  evaluated  include  approaches 
based  on  the  singular  values  of  higher-order  normal  matrices 
for  AR  and  ARMA  models,  and  approaches  based  on  matrix  norms 
applied  to  third-order  cumulants  for  MA  models. 

•  The  performance  of  the  Metford  (1984)  detection  rule  for  non- 
Gaussian  signals  needs  to  be  established.  Alternative 
detection  rules  will  be  identified  and  considered  as 
appropriate.  One  option  is  to  use  non-parametric  detection 
methods  with  the  pseudo- innovations  sequence  (Table  6-1)  . 
Another  option  is  frequency-domain  detection  of  moving 
targets  using  the  bispectrum  (the  Fourier  transform  of  the 
third-order  cumulant  sequence) . 

•  Identification  and  detection  performance  will  be  compared 
with  that  of  other  methods.  This  includes  covariance-based 
methods  using  time  series  as  well  as  state  space  models. 

•  Key  implementation  parameters  will  be  established  for  the 
selected  application ( s ) .  This  includes  the  minimum  required 


122 


channel  output  sequence  duration,  and  the  order  of  the  two 
whitening  filters. 

Real-Time  Processor  Architecture  Design.  A  real-time  architecture 
should  be  developed  for  the  methodology  formulated  in  Phase  I, 
including  identification  of  candidate  hardware.  Such  an 
architecture  should  meet  the  requirements  identified  early  in 
Phase  II  for  the  radar  and  non-defense  applications.  The  very 
large  scale  integration  (VLSI)  systolic  architecture  developed  by 
Liu  and  Jen  (1992)  for  the  IBDA  of  Metford  and  Haykin  (1985) 
provides  a  point  of  departure  since  the  methodology  of  Phase  I  is 
related  to  the  IBDA. 

The  major  computational  task  of  the  methodology  formulated  in 
Phase  I  is  estimation  of  the  third-order  cumulants .  Fortunately, 
the  estimators  for  cumulants  exhibit  a  highly  parallel  structure. 
Stellakis  and  Manolakos  (1993)  have  developed  a  two-stage  VLSI 
architecture  to  estimate  the  third-order  cumulants  of  a  real¬ 
valued  scalar  sequence  which  may  be  extended  to  the  multichannel 
case  (their  architecture  is  applicable  also  to  the  estimation  of 
fourth-order  cumulants) .  In  Phase  II  SSC  will  pursue  extension  of 
the  Stellakis-Manolakos  architecture  to  the  multichannel  case,  and 
explore  alternative  architectures. 
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